TugasSpasial
Kiki Amelia
2026-05-26
#================ 2.1 Input Data Stunting (sesuai file) ================#
library(readxl); library(dplyr); library(stringr)## Warning: package 'dplyr' was built under R version 4.5.2##
## 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)## Warning: package 'ggplot2' was built under R version 4.5.2## 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':
##
## select# 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
TB <- read.csv("DataTBJabar4.csv", header = TRUE, sep = ";", dec = ".")
head(TB)## id Kab.Kota Tahun JumlahPenduduk JumlahKasusTB Ei SMR
## 1 1 Kab Bogor 2020 6067956 16332 13531.512 120.69605
## 2 2 Kab Sukabumi 2020 2526878 5828 5634.925 103.42639
## 3 3 Kab Cianjur 2020 2329635 5068 5195.074 97.55394
## 4 4 Kab Bandung 2020 3773706 7595 8415.346 90.25179
## 5 5 Kab Garut 2020 2644843 3516 5897.987 59.61356
## 6 6 Kab Tasikmalaya 2020 1802165 2620 4018.819 65.19328
## lnSMR AirMinumLayak SanitasiLayak PHBS KepadatanPenduduk HIV
## 1 4.793275 92.13 62.4 52.10 2002 417
## 2 4.638860 79.98 82.9 54.86 657 110
## 3 4.580405 82.73 83.2 64.72 645 189
## 4 4.502603 95.87 75.5 58.90 2050 176
## 5 4.087883 81.86 31.1 60.10 841 5
## 6 4.177356 85.17 91.4 57.20 731 95
## TingkatKemiskinan insidensiTB InsidensiHIV IT IS ST
## 1 7.69 269.1516 6.8721659 1 1 263.3697
## 2 7.09 230.6403 4.3531979 1 2 245.7823
## 3 10.36 217.5448 8.1128589 1 3 242.8362
## 4 6.91 201.2610 4.6638503 1 4 203.8837
## 5 9.98 132.9379 0.1890471 1 5 256.6381
## 6 10.34 145.3807 5.2714374 1 6 220.4697# Join Data Peta dengan Data Stunting
jabar_tb_join <- jabar_sf %>%
left_join(TB, 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(79.951, 159.62, 229.89, 330.03, 469.39, 1171.3)
# compute labels
labels <- c()
brks <- c(79.951, 159.62, 229.89, 330.03, 469.39, 1171.3)
format(brks, scientific=F)## [1] " 79.951" " 159.620" " 229.890" " 330.030" " 469.390" "1171.300"# 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_tb_join$brks <- cut(jabar_tb_join$insidensiTB,
breaks = brks,
include.lowest = TRUE,
labels = labels)
brks_scale <- levels(jabar_tb_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_tb_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("<159,62", "[159.62 – 229.89)", "[229.89 – 330.03)", "[330.03 – 469.39)", ">=469.3"),
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 = "TB (%)") +
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.21] 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) # Untuk as_dgRMatrix_listw jika diperlukan##
## Attaching package: 'spatialreg'## The following objects are masked from 'package:spdep':
##
## get.ClusterOption, get.coresOption, get.mcOption,
## get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOption# 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: NAWB_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)
### Autokorelasi spasial
### Penamaan baris (lebih aman sebelum analisis)
row.names(jabar_sf) <- as.character(seq_len(nrow(jabar_sf)))
### Koordinat titik centroid (pakai sf modern)
CoordK <- st_coordinates(st_centroid(jabar_sf))## Warning: st_centroid assumes attributes are constant over geometries### Plot peta + jaringan tetangga
plot(st_geometry(jabar_sf), axes = TRUE, col = "gray90")
text(CoordK[,1], CoordK[,2],
labels = row.names(jabar_sf),
col = "black", cex = 0.8, pos = 3)
points(CoordK[,1], CoordK[,2],
pch = 19, cex = 0.7, col = "blue")
plot(WB_nb, CoordK, add = TRUE, col = "red", lwd = 2)
##### Open Your Data
DataL<-read.csv("DataTBJabar4.csv", header = TRUE, sep = ";", dec = ".")
head(DataL)## id Kab.Kota Tahun JumlahPenduduk JumlahKasusTB Ei SMR
## 1 1 Kab Bogor 2020 6067956 16332 13531.512 120.69605
## 2 2 Kab Sukabumi 2020 2526878 5828 5634.925 103.42639
## 3 3 Kab Cianjur 2020 2329635 5068 5195.074 97.55394
## 4 4 Kab Bandung 2020 3773706 7595 8415.346 90.25179
## 5 5 Kab Garut 2020 2644843 3516 5897.987 59.61356
## 6 6 Kab Tasikmalaya 2020 1802165 2620 4018.819 65.19328
## lnSMR AirMinumLayak SanitasiLayak PHBS KepadatanPenduduk HIV
## 1 4.793275 92.13 62.4 52.10 2002 417
## 2 4.638860 79.98 82.9 54.86 657 110
## 3 4.580405 82.73 83.2 64.72 645 189
## 4 4.502603 95.87 75.5 58.90 2050 176
## 5 4.087883 81.86 31.1 60.10 841 5
## 6 4.177356 85.17 91.4 57.20 731 95
## TingkatKemiskinan insidensiTB InsidensiHIV IT IS ST
## 1 7.69 269.1516 6.8721659 1 1 263.3697
## 2 7.09 230.6403 4.3531979 1 2 245.7823
## 3 10.36 217.5448 8.1128589 1 3 242.8362
## 4 6.91 201.2610 4.6638503 1 4 203.8837
## 5 9.98 132.9379 0.1890471 1 5 256.6381
## 6 10.34 145.3807 5.2714374 1 6 220.4697DataLY<-data.frame(Case=DataL$insidensiTB)
DataLE<-data.frame(E=DataL$Ei)
DataLX<-data.frame(TingkatKemiskinan=DataL$TingkatKemiskinan, AirMinumLayak=DataL$AirMinumLayak,
SanitasiLayak=DataL$SanitasiLayak, PHBS=DataL$PHBS, KepadatanPenduduk=DataL$KepadatanPenduduk,
InsidensiHIV=DataL$InsidensiHIV)
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 TingkatKemiskinan AirMinumLayak SanitasiLayak PHBS
## 1 269.1516 13531.512 7.69 92.13 62.4 52.10
## 2 230.6403 5634.925 7.09 79.98 82.9 54.86
## 3 217.5448 5195.074 10.36 82.73 83.2 64.72
## 4 201.2610 8415.346 6.91 95.87 75.5 58.90
## 5 132.9379 5897.987 9.98 81.86 31.1 60.10
## 6 145.3807 4018.819 10.34 85.17 91.4 57.20
## KepadatanPenduduk InsidensiHIV
## 1 2002 6.8721659
## 2 657 4.3531979
## 3 645 8.1128589
## 4 2050 4.6638503
## 5 841 0.1890471
## 6 731 5.2714374DATAL1<-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_20=DataL$JumlahKasusTB[1:27]/Ei[1:27]
SMR_21=DataL$JumlahKasusTB[28:54]/Ei[28:54]
SMR_22=DataL$JumlahKasusTB[55:81]/Ei[55:81]
SMR_23=DataL$JumlahKasusTB[82:108]/Ei[82:108]
SMR_24=DataL$JumlahKasusTB[109:135]/Ei[109:135]
SMR_25=DataL$JumlahKasusTB[136:162]/Ei[136:162]
boxplot(SMR_20,SMR_21,SMR_22, SMR_23,SMR_24, SMR_25, main="Boxplot SMR",names=c("2020","2021","2022","2023","2024","2025"))
##### 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 per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.#### 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~TingkatKemiskinan+AirMinumLayak
+SanitasiLayak+PHBS+KepadatanPenduduk+InsidensiHIV+offset(log(E)), family="gaussian", data=FullDataL)
summary(GLM0)##
## Call:
## glm(formula = Case ~ TingkatKemiskinan + AirMinumLayak + SanitasiLayak +
## PHBS + KepadatanPenduduk + InsidensiHIV + offset(log(E)),
## family = "gaussian", data = FullDataL)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 332.262328 232.060388 1.432 0.154217
## TingkatKemiskinan -11.655099 5.789734 -2.013 0.045839 *
## AirMinumLayak 0.988340 2.309453 0.428 0.669278
## SanitasiLayak -2.769743 0.719008 -3.852 0.000171 ***
## PHBS 1.021253 1.164861 0.877 0.381997
## KepadatanPenduduk 0.005016 0.004011 1.251 0.212966
## InsidensiHIV 6.136563 0.762539 8.048 2.06e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 18425.53)
##
## Null deviance: 6223945 on 161 degrees of freedom
## Residual deviance: 2855957 on 155 degrees of freedom
## AIC: 2059.7
##
## Number of Fisher Scoring iterations: 2###(2) Bayesian Poisson Regression
Formula.Cov<-Case~TingkatKemiskinan+AirMinumLayak+SanitasiLayak+PHBS+KepadatanPenduduk+InsidensiHIV
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"))##
## *** 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 valuessummary(Model.Cov)## Time used:
## Pre = 0.908, Running = 0.472, Post = 0.0136, Total = 1.39
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 324.047 227.307 -121.466 324.047 769.560 324.047 0
## TingkatKemiskinan -11.624 5.804 -22.999 -11.624 -0.249 -11.624 0
## AirMinumLayak 1.156 2.277 -3.307 1.156 5.619 1.156 0
## SanitasiLayak -2.765 0.723 -4.181 -2.765 -1.348 -2.765 0
## PHBS 1.033 1.170 -1.259 1.033 3.326 1.033 0
## KepadatanPenduduk 0.005 0.004 -0.003 0.005 0.013 0.005 0
## InsidensiHIV 6.103 0.765 4.603 6.103 7.602 6.103 0
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.00 0.00 0.00 0.00
## 0.975quant mode
## Precision for the Gaussian observations 0.00 0.00
##
## Deviance Information Criterion (DIC) ...............: 2058.03
## Deviance Information Criterion (DIC, saturated) ....: -331.76
## Effective number of parameters .....................: 6.95
##
## Watanabe-Akaike information criterion (WAIC) ...: 2059.91
## Effective number of parameters .................: 8.27
##
## Marginal log-Likelihood: -1093.92
## 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 TingkatKemiskinan AirMinumLayak SanitasiLayak PHBS
## 1 269.1516 13531.512 -0.2262510 -0.3872416 -1.39940415 -1.183567914
## 2 230.6403 5634.925 -0.4466739 -2.5556225 -0.09740126 -0.922716600
## 3 217.5448 5195.074 0.7546310 -2.0648367 -0.07834756 0.009165267
## 4 201.2610 8415.346 -0.5128008 0.2802271 -0.56739255 -0.540890764
## 5 132.9379 5897.987 0.6150298 -2.2201035 -3.38734026 -0.427477150
## 6 145.3807 4018.819 0.7472836 -1.6293758 0.44245359 -0.701560052
## KepadatanPenduduk InsidensiHIV
## 1 -0.4144635 -0.7323041
## 2 -0.7106292 -0.8719440
## 3 -0.7132716 -0.6635259
## 4 -0.4038940 -0.8547228
## 5 -0.6701129 -1.1027850
## 6 -0.6943346 -0.8210411m <- 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)
param.beta = list(prec = list(param = c(1.0e-6, 1.0e-6)))
Formula.Ridge1 = Case ~ f(beta1, TingkatKemiskinan, model="iid", values = c(1:6), hyper = param.beta)+
f(beta2, AirMinumLayak, copy="beta1", fixed=T)+
f(beta3, SanitasiLayak, copy="beta1", fixed=T)+
f(beta4, PHBS, copy="beta1", fixed=T)+
f(beta5, KepadatanPenduduk, copy="beta1", fixed=T)+
f(beta6, InsidensiHIV, 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) 304.6987 15.0427 275.2155 304.6987 334.1819 304.6987 0
## 1 -0.8542 1.7394 -4.2634 -0.8542 2.5550 -0.8542 0
## 2 0.9586 1.7394 -2.4505 0.9586 4.3678 0.9586 0
## 3 -0.3979 1.7393 -3.8069 -0.3979 3.0111 -0.3979 0
## 4 0.3848 1.7394 -3.0242 0.3848 3.7939 0.3848 0
## 5 1.3567 1.7395 -2.0526 1.3567 4.7659 1.3567 0
## 6 1.6555 1.7394 -1.7536 1.6555 5.0647 1.6555 0summary(Model.Ridge1)## Time used:
## Pre = 1.4, Running = 0.487, Post = 0.0187, Total = 1.9
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 304.699 15.043 275.215 304.699 334.182 304.699 0
##
## Random effects:
## Name Model
## beta1 IID model
## beta2 Copy
## beta3 Copy
## beta4 Copy
## beta5 Copy
## beta6 Copy
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.000 0.000 0.000 0.000
## Precision for beta1 0.327 0.017 0.294 0.326
## 0.975quant mode
## Precision for the Gaussian observations 0.000 0.000
## Precision for beta1 0.362 0.325
##
## Deviance Information Criterion (DIC) ...............: 2167.74
## Deviance Information Criterion (DIC, saturated) ....: -201.21
## Effective number of parameters .....................: 1.05
##
## Watanabe-Akaike information criterion (WAIC) ...: 2167.78
## Effective number of parameters .................: 1.09
##
## Marginal log-Likelihood: -1128.52
## 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"))##
## *** 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 valuesCoef.Ridge2 <- Model.Ridge2$summary.random$idx.Z
Coef.Ridge2## ID mean sd 0.025quant 0.5quant 0.975quant
## 1 1 9.685521e-05 0.087195637 -0.1708035 9.685521e-05 0.17099716
## 2 2 5.014630e-03 0.087159747 -0.1658153 5.014630e-03 0.17584459
## 3 3 -1.802395e-02 0.086729746 -0.1880111 -1.802395e-02 0.15196323
## 4 4 9.659250e-03 0.086990440 -0.1608389 9.659250e-03 0.18015738
## 5 5 2.328279e-02 0.002888466 0.0176215 2.328279e-02 0.02894407
## 6 6 5.808095e-02 0.086749476 -0.1119449 5.808095e-02 0.22810680
## mode kld
## 1 9.685521e-05 0.000000e+00
## 2 5.014630e-03 4.467301e-21
## 3 -1.802395e-02 4.349798e-19
## 4 9.659250e-03 8.731996e-20
## 5 2.328279e-02 0.000000e+00
## 6 5.808095e-02 2.306152e-18summary(Model.Ridge2)## Time used:
## Pre = 1.06, Running = 0.399, Post = 0.0114, Total = 1.47
## 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.00 0.00 0.00 0.00
## Precision for idx.Z 139.22 8.38 127.77 137.87
## 0.975quant mode
## Precision for the Gaussian observations 0.00 0.00
## Precision for idx.Z 159.46 131.88
##
## Deviance Information Criterion (DIC) ...............: 2116.53
## Deviance Information Criterion (DIC, saturated) ....: -266.50
## Effective number of parameters .....................: 2.03
##
## Watanabe-Akaike information criterion (WAIC) ...: 2116.99
## Effective number of parameters .................: 2.43
##
## Marginal log-Likelihood: -1099.04
## 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$IT
param.beta = list(prec = list(param = c(1.0e-6, 1.0e-6)))
###(1) Poisson
Formula.Ridge1t = Case ~ f(beta1, TingkatKemiskinan, model="iid", values = c(1:6), hyper = param.beta)+
f(beta2, AirMinumLayak, copy="beta1", fixed=T)+
f(beta3, SanitasiLayak, copy="beta1", fixed=T)+
f(beta4, PHBS, copy="beta1", fixed=T)+
f(beta5, KepadatanPenduduk, copy="beta1", fixed=T)+
f(beta6, InsidensiHIV, 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"))##
## *** 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.Ridge1t <- rbind(Model.Ridge1t$summary.fixed, Model.Ridge1t$summary.random$beta1[,-1])
round(Coeff.Ridge1t,4) ## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 114.8487 29.2532 57.5134 114.8487 172.1839 114.8487 0
## IT2 55.6649 7.5746 40.8189 55.6649 70.5110 55.6649 0
## 1 0.0000 0.0053 -0.0104 0.0000 0.0103 0.0000 0
## 2 0.0000 0.0053 -0.0103 0.0000 0.0104 0.0000 0
## 3 0.0000 0.0053 -0.0104 0.0000 0.0103 0.0000 0
## 4 0.0000 0.0053 -0.0103 0.0000 0.0103 0.0000 0
## 5 0.0000 0.0053 -0.0103 0.0000 0.0104 0.0000 0
## 6 0.0000 0.0053 -0.0103 0.0000 0.0104 0.0000 0summary(Model.Ridge1t)## Time used:
## Pre = 1.56, Running = 0.554, Post = 0.0205, Total = 2.13
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 114.849 29.253 57.513 114.849 172.184 114.849 0
## IT2 55.665 7.575 40.819 55.665 70.511 55.665 0
##
## Random effects:
## Name Model
## beta1 IID model
## IS1 IID model
## beta2 Copy
## beta3 Copy
## beta4 Copy
## beta5 Copy
## beta6 Copy
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.00 0.00 0.00 0.00
## Precision for beta1 33351.91 2815.07 27147.37 33620.44
## Precision for IS1 164607.66 9071.85 149208.33 163902.57
## 0.975quant mode
## Precision for the Gaussian observations 0.00 0.00
## Precision for beta1 37794.70 35795.38
## Precision for IS1 184738.59 161011.21
##
## Deviance Information Criterion (DIC) ...............: 2133.91
## Deviance Information Criterion (DIC, saturated) ....: -246.70
## Effective number of parameters .....................: 1.91
##
## Watanabe-Akaike information criterion (WAIC) ...: 2134.10
## Effective number of parameters .................: 2.07
##
## Marginal log-Likelihood: -1120.63
## 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)')###(1) NonParametric Approach--Temporal Component
### (a) Poisson
param.beta = list(prec = list(param = c(1.0e-6, 1.0e-6)))
Formula.Ridge1tnp = Case ~f(beta1, TingkatKemiskinan, model="iid", values = c(1:6), hyper = param.beta)+
f(beta2, AirMinumLayak, copy="beta1", fixed=T)+
f(beta3, SanitasiLayak, copy="beta1", fixed=T)+
f(beta4, PHBS, copy="beta1", fixed=T)+
f(beta5, KepadatanPenduduk, copy="beta1", fixed=T)+
f(beta6, InsidensiHIV, copy="beta1", fixed=T)+
f(IT1, model = "rw1", hyper = igprior5, cyclic = T, constr = T)
Model.Ridge1tnp <- inla(Formula.Ridge1tnp , 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.Ridge1tnp <- rbind(Model.Ridge1tnp$summary.fixed, Model.Ridge1tnp$summary.random$beta1[,-1])
Coeff.Ridge1tnp## mean sd 0.025quant 0.5quant 0.975quant mode
## (Intercept) 311.753168 0.01062936 311.732335 311.753168 311.774001 311.753168
## 1 -8.481061 0.01596886 -8.512359 -8.481061 -8.449762 -8.481061
## 2 -6.799885 0.01297867 -6.825322 -6.799885 -6.774447 -6.799885
## 3 -27.379211 0.01193710 -27.402607 -27.379211 -27.355815 -27.379211
## 4 -2.067205 0.01240205 -2.091512 -2.067205 -2.042897 -2.067205
## 5 53.618700 0.01854058 53.582361 53.618700 53.655039 53.618700
## 6 94.102701 0.01407032 94.075123 94.102701 94.130278 94.102701
## kld
## (Intercept) 0.000000e+00
## 1 4.050140e-13
## 2 0.000000e+00
## 3 2.993357e-10
## 4 2.397214e-12
## 5 0.000000e+00
## 6 3.765723e-09summary(Model.Ridge1tnp)## Time used:
## Pre = 1.77, Running = 0.609, Post = 0.0307, Total = 2.41
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 311.753 0.011 311.732 311.753 311.774 311.753 0
##
## Random effects:
## Name Model
## beta1 IID model
## IT1 RW1 model
## beta2 Copy
## beta3 Copy
## beta4 Copy
## beta5 Copy
## beta6 Copy
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 54.63 0.00 54.63 54.63
## Precision for beta1 54.60 0.00 54.60 54.60
## Precision for IT1 54.65 0.00 54.65 54.65
## 0.975quant mode
## Precision for the Gaussian observations 54.63 54.63
## Precision for beta1 54.60 54.60
## Precision for IT1 54.65 54.65
##
## Deviance Information Criterion (DIC) ...............: -92800457.84
## Deviance Information Criterion (DIC, saturated) ....: -92800451.31
## Effective number of parameters .....................: -92800672.63
##
## Watanabe-Akaike information criterion (WAIC) ...: 156208.62
## Effective number of parameters .................: 72333.33
##
## Marginal log-Likelihood: -48205499.77
## 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)')### (3) Spatial component---Besag
### (a) Poisson
Formula.Ridge1s = Case ~f(beta1, TingkatKemiskinan, model="iid", values = c(1:6), hyper = param.beta)+
f(beta2, AirMinumLayak, copy="beta1", fixed=T)+
f(beta3, SanitasiLayak, copy="beta1", fixed=T)+
f(beta4, PHBS, copy="beta1", fixed=T)+
f(beta5, KepadatanPenduduk, copy="beta1", fixed=T)+
f(beta6, InsidensiHIV, 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) ## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 311.7532 0.0127 311.7282 311.7532 311.7781 311.7532 0
## 1 -138.0259 0.0424 -138.1089 -138.0259 -137.9429 -138.0259 0
## 2 18.5679 0.0276 18.5139 18.5679 18.6220 18.5679 0
## 3 -23.1115 0.0172 -23.1453 -23.1115 -23.0778 -23.1115 0
## 4 19.6099 0.0219 19.5670 19.6099 19.6528 19.6099 0
## 5 7.5256 0.0455 7.4365 7.5256 7.6148 7.5256 0
## 6 42.0325 0.0252 41.9831 42.0325 42.0819 42.0325 0summary(Model.Ridge1s )## Time used:
## Pre = 1.91, Running = 0.54, Post = 0.058, Total = 2.5
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 311.753 0.013 311.728 311.753 311.778 311.753 0
##
## Random effects:
## Name Model
## beta1 IID model
## IS1 Besags ICAR model
## beta2 Copy
## beta3 Copy
## beta4 Copy
## beta5 Copy
## beta6 Copy
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 38.08 0.00 38.08 38.08
## Precision for beta1 105.50 0.00 105.50 105.50
## Precision for IS1 10.69 0.00 10.69 10.69
## 0.975quant mode
## Precision for the Gaussian observations 38.08 38.08
## Precision for beta1 105.50 105.50
## Precision for IS1 10.69 10.69
##
## Deviance Information Criterion (DIC) ...............: -67454905.05
## Deviance Information Criterion (DIC, saturated) ....: -67454896.20
## Effective number of parameters .....................: -67455015.92
##
## Watanabe-Akaike information criterion (WAIC) ...: 71605.98
## Effective number of parameters .................: 30094.39
##
## Marginal log-Likelihood: -38673404.64
## 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)')### (4) Spatial component---BYM
## (a) Poisson
Formula.Ridge1s2 = Case ~ f(beta1, TingkatKemiskinan, model="iid", values = c(1:6), hyper = param.beta)+
f(beta2, AirMinumLayak, copy="beta1", fixed=T)+
f(beta3, SanitasiLayak, copy="beta1", fixed=T)+
f(beta4, PHBS, copy="beta1", fixed=T)+
f(beta5, KepadatanPenduduk, copy="beta1", fixed=T)+
f(beta6, InsidensiHIV, 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) ## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 311.7532 0.0106 311.7323 311.7532 311.7740 311.7532 0
## 1 -156.2423 0.0389 -156.3184 -156.2423 -156.1661 -156.2423 0
## 2 22.2273 0.0242 22.1798 22.2273 22.2748 22.2273 0
## 3 -22.8033 0.0147 -22.8320 -22.8033 -22.7746 -22.8033 0
## 4 16.5325 0.0196 16.4941 16.5325 16.5708 16.5325 0
## 5 -13.6366 0.0418 -13.7186 -13.6366 -13.5546 -13.6366 0
## 6 48.1958 0.0213 48.1540 48.1958 48.2376 48.1958 0summary(Model.Ridge1s2)## Time used:
## Pre = 1.89, Running = 0.606, Post = 0.0235, Total = 2.52
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 311.753 0.011 311.732 311.753 311.774 311.753 0
##
## Random effects:
## Name Model
## beta1 IID model
## IS1 Besags ICAR model
## IS2 IID model
## beta2 Copy
## beta3 Copy
## beta4 Copy
## beta5 Copy
## beta6 Copy
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 54.59 0.00 54.59 54.59
## Precision for beta1 54.61 0.00 54.61 54.61
## Precision for IS1 54.60 0.00 54.60 54.60
## Precision for IS2 54.62 0.00 54.62 54.62
## 0.975quant mode
## Precision for the Gaussian observations 54.59 54.59
## Precision for beta1 54.61 54.61
## Precision for IS1 54.60 54.60
## Precision for IS2 54.62 54.62
##
## Deviance Information Criterion (DIC) ...............: -93884858.19
## Deviance Information Criterion (DIC, saturated) ....: -93884851.71
## Effective number of parameters .....................: -93884963.23
##
## Watanabe-Akaike information criterion (WAIC) ...: 137235.12
## Effective number of parameters .................: 62842.01
##
## Marginal log-Likelihood: -54771146.27
## 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)')### (5) Spatial + Temporal component---Besag Parametric
Formula.Ridge1nbst = Case ~ f(beta1, TingkatKemiskinan, model="iid", values = c(1:6), hyper = param.beta)+
f(beta2, AirMinumLayak, copy="beta1", fixed=T)+
f(beta3, SanitasiLayak, copy="beta1", fixed=T)+
f(beta4, PHBS, copy="beta1", fixed=T)+
f(beta5, KepadatanPenduduk, copy="beta1", fixed=T)+
f(beta6, InsidensiHIV, 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"))##
## *** 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.Ridge1nbst <- rbind(Model.Ridge1nbst$summary.fixed, Model.Ridge1nbst$summary.random$beta1[,-1])
round(Coeff.Ridge1nbst,4) ## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 114.8605 29.3131 57.4079 114.8605 172.3132 114.8605 0
## IT2 55.6588 7.5905 40.7817 55.6588 70.5359 55.6588 0
## 1 -0.0411 0.3716 -0.7695 -0.0411 0.6873 -0.0411 0
## 2 0.0470 0.3716 -0.6814 0.0470 0.7754 0.0470 0
## 3 -0.0218 0.3716 -0.7502 -0.0218 0.7065 -0.0218 0
## 4 0.0067 0.3716 -0.7217 0.0067 0.7351 0.0067 0
## 5 0.0799 0.3716 -0.6484 0.0799 0.8083 0.0799 0
## 6 0.0865 0.3716 -0.6419 0.0865 0.8148 0.0865 0summary(Model.Ridge1nbst)## Time used:
## Pre = 1.79, Running = 0.705, Post = 0.0273, Total = 2.52
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 114.861 29.313 57.408 114.861 172.313 114.861 0
## IT2 55.659 7.590 40.782 55.659 70.536 55.659 0
##
## Random effects:
## Name Model
## beta1 IID model
## IS1 Besags ICAR model
## IS3 IID model
## beta2 Copy
## beta3 Copy
## beta4 Copy
## beta5 Copy
## beta6 Copy
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.00 0.000 0.00 0.00
## Precision for beta1 7.50 0.407 6.87 7.45
## Precision for IS1 44.81 2.850 40.43 44.45
## Precision for IS3 51.92 2.090 48.27 51.79
## 0.975quant mode
## Precision for the Gaussian observations 0.00 0.00
## Precision for beta1 8.45 7.24
## Precision for IS1 51.44 43.02
## Precision for IS3 56.48 51.22
##
## Deviance Information Criterion (DIC) ...............: 2133.62
## Deviance Information Criterion (DIC, saturated) ....: -246.80
## Effective number of parameters .....................: 1.91
##
## Watanabe-Akaike information criterion (WAIC) ...: 2133.80
## Effective number of parameters .................: 2.07
##
## Marginal log-Likelihood: -1139.18
## 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)')### (6) Spatial + Temporal component---Besag NonParametric
## (a) Poisson
Formula.Ridge1stnp = Case ~ f(beta1, TingkatKemiskinan, model="iid", values = c(1:6), hyper = param.beta)+
f(beta2, AirMinumLayak, copy="beta1", fixed=T)+
f(beta3, SanitasiLayak, copy="beta1", fixed=T)+
f(beta4, PHBS, copy="beta1", fixed=T)+
f(beta5, KepadatanPenduduk, copy="beta1", fixed=T)+
f(beta6, InsidensiHIV, copy="beta1", fixed=T)+
f(IS1,model = "besag", graph = W1, constr = T, hyper = hcprior)+
f(IT1, model = "rw1", hyper = igprior5, constr = T, cyclic=T)
Model.Ridge1stnp <- inla(Formula.Ridge1stnp , 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.Ridge1stnp <- rbind(Model.Ridge1stnp$summary.fixed, Model.Ridge1stnp$summary.random$beta1[,-1])
round(Coeff.Ridge1stnp,4) ## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 304.7926 14.9423 275.5062 304.7926 334.0790 304.7926 0
## 1 -0.0021 0.0848 -0.1684 -0.0021 0.1642 -0.0021 0
## 2 0.0024 0.0848 -0.1639 0.0024 0.1687 0.0024 0
## 3 -0.0010 0.0848 -0.1672 -0.0010 0.1653 -0.0010 0
## 4 0.0009 0.0848 -0.1653 0.0009 0.1672 0.0009 0
## 5 0.0033 0.0848 -0.1630 0.0033 0.1696 0.0033 0
## 6 0.0040 0.0848 -0.1623 0.0040 0.1703 0.0040 0summary(Model.Ridge1stnp)## Time used:
## Pre = 1.56, Running = 0.827, Post = 0.0283, Total = 2.42
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 304.793 14.942 275.506 304.793 334.079 304.793 0
##
## Random effects:
## Name Model
## beta1 IID model
## IS1 Besags ICAR model
## IT1 RW1 model
## beta2 Copy
## beta3 Copy
## beta4 Copy
## beta5 Copy
## beta6 Copy
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.000 0.000 0.000 0.000
## Precision for beta1 138.990 5.880 127.750 138.874
## Precision for IS1 107.697 5.378 97.580 107.534
## Precision for IT1 0.367 0.018 0.337 0.366
## 0.975quant mode
## Precision for the Gaussian observations 0.000 0.000
## Precision for beta1 150.898 138.661
## Precision for IS1 118.753 107.148
## Precision for IT1 0.406 0.361
##
## Deviance Information Criterion (DIC) ...............: 2171.80
## Deviance Information Criterion (DIC, saturated) ....: -197.84
## Effective number of parameters .....................: 0.984
##
## Watanabe-Akaike information criterion (WAIC) ...: 2171.84
## Effective number of parameters .................: 1.02
##
## Marginal log-Likelihood: -1168.80
## 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)')### (7) Spatial + Temporal component---Interaction
Formula.Ridge1nbstInt = Case ~
f(beta1, TingkatKemiskinan, model="iid", values=c(1:6), hyper=param.beta)+
f(beta2, AirMinumLayak, copy="beta1", fixed=TRUE)+
f(beta3, SanitasiLayak, copy="beta1", fixed=TRUE)+
f(beta4, PHBS, copy="beta1", fixed=TRUE)+
f(beta5, KepadatanPenduduk, copy="beta1", fixed=TRUE)+
f(beta6, InsidensiHIV, copy="beta1", fixed=TRUE)+
f(IS1,
model="besag",
graph=W1,
constr=TRUE,
scale.model=TRUE,
hyper=hcprior,
group=IT2,
control.group=list(
model="rw1",
cyclic=FALSE,
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) ## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 311.7532 0.0152 311.7233 311.7532 311.7830 311.7532 0
## 1 -384.6991 0.0956 -384.8866 -384.6991 -384.5117 -384.6991 0
## 2 24.6238 0.0534 24.5191 24.6238 24.7286 24.6238 0
## 3 -1.1835 0.0291 -1.2405 -1.1835 -1.1265 -1.1835 0
## 4 28.6055 0.0390 28.5291 28.6055 28.6819 28.6055 0
## 5 33.1653 0.1763 32.8198 33.1653 33.5108 33.1653 0
## 6 40.9516 0.0485 40.8566 40.9516 41.0466 40.9516 0summary(Model.Ridge1nbstInt)## Time used:
## Pre = 1.77, Running = 0.608, Post = 0.0241, Total = 2.4
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 311.753 0.015 311.723 311.753 311.783 311.753 0
##
## Random effects:
## Name Model
## beta1 IID model
## IS1 Besags ICAR model
## beta2 Copy
## beta3 Copy
## beta4 Copy
## beta5 Copy
## beta6 Copy
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 26.57 0.00 26.57 26.57
## Precision for beta1 30.81 0.00 30.81 30.81
## Precision for IS1 36.91 0.00 36.91 36.91
## 0.975quant mode
## Precision for the Gaussian observations 26.57 26.57
## Precision for beta1 30.81 30.81
## Precision for IS1 36.91 36.91
##
## Deviance Information Criterion (DIC) ...............: -10307052.21
## Deviance Information Criterion (DIC, saturated) ....: -10307047.07
## Effective number of parameters .....................: -10307105.84
##
## Watanabe-Akaike information criterion (WAIC) ...: 12911.24
## Effective number of parameters .................: 754.74
##
## Marginal log-Likelihood: -9764712.84
## 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)')ST_20=Model.Ridge1nbstInt$summary.fitted.values$mean[1:27]
ST_21=Model.Ridge1nbstInt$summary.fitted.values$mean[28:54]
ST_22=Model.Ridge1nbstInt$summary.fitted.values$mean[55:81]
ST_23=Model.Ridge1nbstInt$summary.fitted.values$mean[82:108]
ST_24=Model.Ridge1nbstInt$summary.fitted.values$mean[109:135]
ST_25=Model.Ridge1nbstInt$summary.fitted.values$mean[136:162]
boxplot(ST_20,ST_21,ST_22, ST_23, ST_24,ST_25, main="Boxplot ST",names=c("2020","2021","2022","2023","2024","2025"))
ST=c(ST_20,ST_21,ST_22, ST_23, ST_24,ST_25)
TB <- read.csv("DataTBJabar4.csv", header = TRUE, sep = ";", dec = ".")
head(TB)## id Kab.Kota Tahun JumlahPenduduk JumlahKasusTB Ei SMR
## 1 1 Kab Bogor 2020 6067956 16332 13531.512 120.69605
## 2 2 Kab Sukabumi 2020 2526878 5828 5634.925 103.42639
## 3 3 Kab Cianjur 2020 2329635 5068 5195.074 97.55394
## 4 4 Kab Bandung 2020 3773706 7595 8415.346 90.25179
## 5 5 Kab Garut 2020 2644843 3516 5897.987 59.61356
## 6 6 Kab Tasikmalaya 2020 1802165 2620 4018.819 65.19328
## lnSMR AirMinumLayak SanitasiLayak PHBS KepadatanPenduduk HIV
## 1 4.793275 92.13 62.4 52.10 2002 417
## 2 4.638860 79.98 82.9 54.86 657 110
## 3 4.580405 82.73 83.2 64.72 645 189
## 4 4.502603 95.87 75.5 58.90 2050 176
## 5 4.087883 81.86 31.1 60.10 841 5
## 6 4.177356 85.17 91.4 57.20 731 95
## TingkatKemiskinan insidensiTB InsidensiHIV IT IS ST
## 1 7.69 269.1516 6.8721659 1 1 263.3697
## 2 7.09 230.6403 4.3531979 1 2 245.7823
## 3 10.36 217.5448 8.1128589 1 3 242.8362
## 4 6.91 201.2610 4.6638503 1 4 203.8837
## 5 9.98 132.9379 0.1890471 1 5 256.6381
## 6 10.34 145.3807 5.2714374 1 6 220.4697# Join Data Peta dengan Data Stunting
jabar_tb_join <- jabar_sf %>%
left_join(TB, 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(79.951, 159.62, 229.89, 330.03, 469.39, 1171.3)
# compute labels
labels <- c()
brks <- c(79.951, 159.62, 229.89, 330.03, 469.39, 1171.3)
format(brks, scientific=F)## [1] " 79.951" " 159.620" " 229.890" " 330.030" " 469.390" "1171.300"# 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_tb_join$brks <- cut(jabar_tb_join$ST,
breaks = brks,
include.lowest = TRUE,
labels = labels)
brks_scale <- levels(jabar_tb_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_tb_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("<159,62", "[159.62 – 229.89)", "[229.89 – 330.03)", "[330.03 – 469.39)", ">=469.3"),
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 = "TB (%)") +
theme(legend.position="right", text = element_text(size=25)) +
theme(text = element_text(size=13))
A
### (8) Spatial + Temporal component---Interaction Poisson
Formula.Ridge1stInt = Case ~
f(beta1, TingkatKemiskinan, model="iid", values=c(1:6), hyper=param.beta)+
f(beta2, AirMinumLayak, copy="beta1", fixed=TRUE)+
f(beta3, SanitasiLayak, copy="beta1", fixed=TRUE)+
f(beta4, PHBS, copy="beta1", fixed=TRUE)+
f(beta5, KepadatanPenduduk, copy="beta1", fixed=TRUE)+
f(beta6, InsidensiHIV, copy="beta1", fixed=TRUE)+
f(IS1,
model="besag",
graph=W1,
constr=TRUE,
hyper=hcprior,
group=IT2,
control.group=list(
model="rw1",
cyclic=FALSE,
hyper=igprior5
)
)
Model.Ridge1stInt <- inla(Formula.Ridge1stInt , 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.Ridge1stInt <- rbind(Model.Ridge1stInt$summary.fixed, Model.Ridge1stInt$summary.random$beta1[,-1])
round(Coeff.Ridge1stInt,4) ## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 308.6728 9.9402 289.1904 308.6728 328.1552 308.6728 0
## 1 -39.6544 18.0993 -75.1285 -39.6544 -4.1804 -39.6544 0
## 2 11.2983 14.5193 -17.1589 11.2983 39.7556 11.2983 0
## 3 -40.7401 11.5748 -63.4264 -40.7401 -18.0538 -40.7401 0
## 4 13.1157 13.3213 -12.9936 13.1157 39.2249 13.1157 0
## 5 20.5195 20.7514 -20.1526 20.5195 61.1915 20.5195 0
## 6 99.6142 14.8059 70.5951 99.6142 128.6333 99.6142 0summary(Model.Ridge1stInt) ## Time used:
## Pre = 1.5, Running = 1.68, Post = 0.0289, Total = 3.21
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 308.673 9.94 289.19 308.673 328.155 308.673 0
##
## Random effects:
## Name Model
## beta1 IID model
## IS1 Besags ICAR model
## beta2 Copy
## beta3 Copy
## beta4 Copy
## beta5 Copy
## beta6 Copy
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.000 0.000 0.000 0.000
## Precision for beta1 0.000 0.000 0.000 0.000
## Precision for IS1 0.234 0.013 0.214 0.232
## 0.975quant mode
## Precision for the Gaussian observations 0.000 0.000
## Precision for beta1 0.000 0.000
## Precision for IS1 0.264 0.226
##
## Deviance Information Criterion (DIC) ...............: 2052.43
## Deviance Information Criterion (DIC, saturated) ....: -333.78
## Effective number of parameters .....................: 15.38
##
## Watanabe-Akaike information criterion (WAIC) ...: 2054.29
## Effective number of parameters .................: 15.78
##
## Marginal log-Likelihood: -1032.21
## 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)')ST_20=Model.Ridge1stInt$summary.fitted.values$mean[1:27]
ST_21=Model.Ridge1stInt$summary.fitted.values$mean[28:54]
ST_22=Model.Ridge1stInt$summary.fitted.values$mean[55:81]
ST_23=Model.Ridge1stInt$summary.fitted.values$mean[82:108]
ST_24=Model.Ridge1stInt$summary.fitted.values$mean[109:135]
ST_25=Model.Ridge1stInt$summary.fitted.values$mean[136:162]
boxplot(ST_20,ST_21,ST_22, ST_23, ST_24,ST_25, main="Boxplot ST",names=c("2020","2021","2022","2023","2024","2025"))
ST=c(ST_20,ST_21,ST_22, ST_23, ST_24,ST_25)(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)