Preparation for Master Thesis
Dila
2024-12-24
#set working directory
setwd("/Users/dila.azuri/Documents/Magister/Tesis/SPDE/Syntax")
#====================== 1. Install and read the Library ==================================
#install.packages("INLA", repos = "https://inla.r-inla-download.org/R/stable", dep = TRUE)
#install.packages("devtools")
#devtools::install_github(repo = "https://github.com/hrue/r-inla", ref = "stable", subdir = "rinla", build = FALSE)
library("devtools")
library(ggplot2)
library(sf)
library(gstat)
library(sp)
library(INLA) #INLA
library(writexl) #save to xls
library(dplyr) #for pivot in statistics descriptive
library(plotly) # for surface plot
library(akima) #for smoothness interpolation result
library(cowplot) #gabungin plot
library(patchwork) #label title di ggplot
#================================== 2. Input Data ==================================
#====================== 2.1 Earthquake data
Earthquakes<-read.csv("INLA - Gempa.csv",sep=";")
head(Earthquakes)
#======= Data 2020
Earthquakes20=Earthquakes[Earthquakes$Tahun==2020,]
head(Earthquakes20)
str(Earthquakes20)
#======= Data Preprocessing
Earthquakes20$LAT <- gsub(",", ".", Earthquakes20$LAT)
Earthquakes20$LAT <- as.numeric(Earthquakes20$LAT)
Earthquakes20$LON <- gsub(",", ".", Earthquakes20$LON)
Earthquakes20$LON <- as.numeric(Earthquakes20$LON)
Earthquakes20$MAG <- gsub(",", ".", Earthquakes20$MAG)
Earthquakes20$MAG <- as.numeric(Earthquakes20$MAG)
colnames(Earthquakes20)[colnames(Earthquakes20) == "Bulan"] <- "Month"
str(Earthquakes20)
head(Earthquakes20)
#====================== 2.2 Sumatra Map
Indonesia<-readRDS('gadm36_IDN_1_sp.rds')
Indonesia$NAME_1
Sumatra<-Indonesia[Indonesia$NAME_1 %in%c("Kepulauan Riau","Sumatera Selatan","Lampung","Sumatera Utara",
"Bangka Belitung","Jambi","Riau","Bengkulu","Sumatera Barat","Aceh"),]
Sumatra_sf <- st_as_sf(Sumatra)
#====================== 2.3 Fault Map:5000.shp
patahan = st_read(file.choose())
sumatera_patahan <- st_bbox(c(xmin = 95, ymin = -6, xmax = 106, ymax = 6), crs = st_crs(patahan))
patahan_sumatera <- st_crop(patahan, sumatera_patahan)
#====================== 2.4 Visualization
Month_names <- c( "1"="January",
"2"="February",
"3"="March",
"4"="April",
"5"="May",
"6"="June",
"7"="July",
"8"="August",
"9"="September",
"10"="October",
"11"="November",
"12"="December")
#================================== 3. New Data set ==================================
# Create a unique list of all Longitude and Latitude combinations
locations <- unique(Earthquakes20[, c("LON", "LAT")])
# Create a grid of all combinations of Longitude, Latitude, and Time
all_combinations <- merge(locations, data.frame(Month = unique(Earthquakes20$Month)))
head(all_combinations)
# Merge with original Earthquakes20
final_Earthquakes20 <- merge(all_combinations, Earthquakes20,
by = c("LON", "LAT", "Month"), all.x = TRUE)
head(final_Earthquakes20)
# Fill NA for Locations Missing Data at Specific Times
final_Earthquakes20$MAG[is.na(final_Earthquakes20$MAG)] <- NA
final_Earthquakes20$DEPTH[is.na(final_Earthquakes20$DEPTH)] <- NA
final_Earthquakes20$Year[is.na(final_Earthquakes20$Year)] <- 2020
# Construct final dataframe
Earthquakes20<-data.frame(LON=final_Earthquakes20$LON,LAT=final_Earthquakes20$LAT,
Month=final_Earthquakes20$Month, MAG=final_Earthquakes20$MAG,
DEPTH=final_Earthquakes20$DEPTH)
head(Earthquakes20)
#================================== 4. CREATE GRID ==================================
# Geographic coordinate data (longitude and latitude)
polygon_coords <- data.frame(
longitude = c(94, 94, 98, 108, 102, 94),
latitude = c(2, 6, 6, -7, -7, 2)
)
polygon_sf <- st_as_sf(polygon_coords, coords = c("longitude", "latitude"), crs = 4326)
# Convert to UTM (e.g., UTM Zone 48N)
polygon_utm <- st_transform(polygon_sf, crs = 32648)
utm_coords <- as.data.frame(st_coordinates(polygon_utm))
# Membuat data frame untuk utm_coords jika belum ada
utm_coords_df <- data.frame(longitude = utm_coords$X, latitude = utm_coords$Y)
#====================== 4.1 GRID: 10km x 10km
# Define the interval between grids. Update if needed
interval1 <- 10000
grid_points1 <- expand.grid(
longitude = seq(min(utm_coords$X), max(utm_coords$X), by = interval1),
latitude = seq(min(utm_coords$Y), max(utm_coords$Y), by = interval1)
)
inside1 <- point.in.polygon(grid_points1$longitude, grid_points1$latitude,
utm_coords$X, utm_coords$Y)
inside_points1 <- grid_points1[inside1 == 1, ]
## Square grid with constant interval
coordinates(inside_points1)<-~longitude+latitude
gridded(inside_points1) = TRUE
coordsGrid.UTM1<-inside_points1 ## Get the UTM grid location
UTM<-CRS("+proj=utm +zone=48 ellps=WGS84")
Sumatra_UTM <- spTransform(Sumatra, UTM)
coordsGrid.UTM_df1 <- as.data.frame(coordsGrid.UTM1)
Sumatra_UTM_sf <- st_as_sf(Sumatra_UTM)
#============================ 5. Prepare the data in the UTM system ===========================
st_crs(Sumatra_sf)
Sumatra_sf_utm1 <- st_transform(Sumatra_sf, crs = 32648)
## Observed
Observed<-data.frame(X=Earthquakes20$LON,Y=Earthquakes20$LAT)
coordinates(Observed) <- ~X+Y # Asumsikan kolom X dan Y adalah kolom lon dan lat
proj4string(Observed) <- CRS("+proj=longlat +datum=WGS84") # Set CRS ke WGS84 (EPSG: 4326)
Observed_utm <- spTransform(Observed, CRS("+proj=utm +zone=48 +datum=WGS84 +units=m +no_defs"))
Observed_utm1<-coordinates(Observed_utm)
colnames(Observed_utm1)<-c("X","Y")
Earthquakes20$UTM_x<-Observed_utm1[,1]
Earthquakes20$UTM_y<-Observed_utm1[,2]
Earthquakes20$IT<-Earthquakes20$Month
head(Earthquakes20)1. Gaussian
#cek distribusi gaussian engga
#1. histogram
hist(Earthquakes20$MAG, probability = TRUE, main = "Histogram dengan Kurva Normal")
lines(density(Earthquakes20$MAG), col = "blue", lwd = 2) # Kernel density
curve(dnorm(x, mean = mean(Earthquakes20$MAG), sd = sd(Earthquakes20$MAG)), add = TRUE, col = "red", lwd = 2) # Kurva normal
#2. qqplot
qqnorm(Earthquakes20$MAG, main = "Q-Q Plot")
qqline(Earthquakes20$MAG, col = "red", lwd = 2) # Garis referensi
#3. Kolmogorov smirnov: terima Ho artinya normal
ks.test(Earthquakes20$MAG, "pnorm", mean = mean(Earthquakes20$MAG), sd = sd(Earthquakes20$MAG))
#
# Asymptotic one-sample Kolmogorov-Smirnov test
#
#data: Earthquakes20$MAG
#D = 0.095553, p-value = 0.2726
#alternative hypothesis: two-sided2. Mesh
par(mfrow = c(2, 3))
#Mesh 0
Earthquakes_mesh0 <- inla.mesh.2d(
loc = cbind(Earthquakes20$UTM_x, Earthquakes20$UTM_y), # UTM coordinates of earthquakes
max.edge = 100000
)
num_vertices0 <- Earthquakes_mesh0$n
plot(Earthquakes_mesh0, col = "lightgray", main = "Mesh 0: Max Edge = 80 km")
points(Observed_utm1,col="red")
# Mesh 1
Earthquakes_mesh1 <- inla.mesh.2d(
loc = cbind(Earthquakes20$UTM_x, Earthquakes20$UTM_y), # UTM coordinates of earthquakes
max.edge = c(100000, 1000000) # Maximum triangle edge length
)
num_vertices1 <- Earthquakes_mesh1$n
plot(Earthquakes_mesh1, col = "lightgray", main = "Mesh 1: Max Edge = 80 km to 1M km")
points(Observed_utm1,col="red")
# Mesh 2
Earthquakes_mesh2 <- inla.mesh.2d(
loc = cbind(Earthquakes20$UTM_x, Earthquakes20$UTM_y), # UTM coordinates of earthquakes
max.edge = c(100000, 1300000) # Maximum triangle edge length
)
num_vertices2 <- Earthquakes_mesh2$n
plot(Earthquakes_mesh2, col = "lightgray", main = "Mesh 2: Max Edge = 80 km to 1.3M km")
points(Observed_utm1,col="red")
# Mesh 3
Earthquakes_mesh3 <- inla.mesh.2d(
loc = cbind(Earthquakes20$UTM_x, Earthquakes20$UTM_y), # UTM coordinates of earthquakes
offset = c(50000, 200000), # Min and max distance to mesh boundary
max.edge = c(100000, 1000000) # Maximum triangle edge length
)
num_vertices3 <- Earthquakes_mesh3$n
plot(Earthquakes_mesh3, col = "lightgray", main = "Mesh 3: Offset = 50 km to 100 km")
points(Observed_utm1,col="red")
# Mesh 4
Earthquakes_mesh4 <- inla.mesh.2d(
loc = cbind(Earthquakes20$UTM_x, Earthquakes20$UTM_y), # UTM coordinates of earthquakes
offset = c(200000, 50000), # Min and max distance to mesh boundary
max.edge = c(100000, 1000000) # Maximum triangle edge length
)
num_vertices4 <- Earthquakes_mesh4$n
plot(Earthquakes_mesh4, col = "lightgray", main = "Mesh 4: Offset = 100 km to 50 km")
points(Observed_utm1,col="red")Mesh dari kiri ke kanan dan bawah. Cek Buku Blangiardo 2015 hal 201. Pakenya Mesh3
3. Modeling
coordinates.allyear<-as.matrix(data.frame(x=Earthquakes20$UTM_x,y=Earthquakes20$UTM_y))
control <- list(
results = list(return.marginals.random = TRUE, return.marginals.predictor=TRUE),
compute = list(hyperpar=TRUE, return.marginals.predictor=TRUE, return.marginals=TRUE, dic=TRUE, mlik = TRUE, cpo = TRUE,
po = TRUE, waic=TRUE))
#=========================== Grid: 10km * 10km
Grid.CoordUTM1<-as.data.frame(coordsGrid.UTM1)
m1<-nrow(Grid.CoordUTM1)
Grid.CoordUTM1<-as.matrix(Grid.CoordUTM1)
GroupTime1<-rep(c(1:12),each=m1)
All.Grid1<-kronecker(matrix(1, 12, 1), Grid.CoordUTM1)# 1. Create the mesh (assuming mesh has been defined as Earthquakes_mesh)
A_est <- inla.spde.make.A(mesh=Earthquakes_mesh3,loc=coordinates.allyear,
group=Earthquakes20$Month, n.group=12)
dim(A_est)
# 2. Define the SPDE model (using default priors, no specific priors set)
Earthquakes_spde <- inla.spde2.pcmatern(
mesh = Earthquakes_mesh3,
prior.range = c(1000, 0.5), # default
prior.sigma = c(1, 0.5) # default
)
# 3. Create index for spatial field
s_index <- inla.spde.make.index(name = "spatial.field",
n.spde = Earthquakes_spde$n.spde, n.group = 12 # Example for grouping time or space as needed
)
# 4. Create stack for estimation (training data)
stack_est <- inla.stack(
data = list(MAG = log(Earthquakes20$MAG)), # Example log-transformed MAG data
A = list(1, A_est), # Adjust A_est as required for your model
effects = list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time = Earthquakes20$Month, UTM_x = Earthquakes20$UTM_x / 1000,
UTM_y = Earthquakes20$UTM_y / 1000, DEPTH = Earthquakes20$DEPTH),
s_index), tag = "est")
# 5. Create stack for prediction (unobserved or new data)
A_pred <- inla.spde.make.A(mesh=Earthquakes_mesh3,loc=as.matrix(All.Grid1),
group=GroupTime1, #selected day for prediction
n.group=12)
Cov<-data.frame(UTM_x=All.Grid1[,1]/1000,UTM_y=All.Grid1[,2]/1000,DEPTH=NA,Time=GroupTime1)
stack_pred <- inla.stack(
data = list(MAG = NA), A = list(1, A_pred), # Adjust A_pred5 as required for prediction
effects = list(data.frame(
Intercept = rep(1, nrow(All.Grid1)),Time = GroupTime1, # Assuming GroupTime5 is defined
UTM_x = Cov$UTM_x,UTM_y = Cov$UTM_y,DEPTH = Cov$DEPTH), s_index # Spatial index for SPDE
), tag = "pred")
# Combine both estimation and prediction stacks
stack <- inla.stack(stack_est, stack_pred) 3.1 Spatiotemporal: Spatial with time dengan interaksi RW 1
# 6. Define the formula (specifying fixed effects and SPDE random field)
# hal 22: https://sites.stat.washington.edu/peter/591/Lindgren.pdf
formula <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model = Earthquakes_spde, group = spatial.field.group,
control.group = list(model = "rw1"))
# 7. Run the INLA model (without any specific priors defined)
output <- inla(formula,data=inla.stack.data(stack, spde=Earthquakes_spde),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack), compute=TRUE))
# 8. View summary of the output
summary(output)Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data,
quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials =
Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
verbose, ", " lincomb = lincomb, selection = selection, control.compute = control.compute, ", "
control.predictor = control.predictor, control.family = control.family, ", " control.inla =
control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert =
control.expert, ", " control.hazard = control.hazard, control.lincomb = control.lincomb, ", "
control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso =
control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg =
inla.arg, num.threads = num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug = debug, .parent.frame =
.parent.frame)" )
Time used:
Pre = 2.23, Running = 91.4, Post = 5.31, Total = 98.9
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
Intercept 0.376 23.307 -48.558 0.495 49.020 0.459 0
UTM_x -0.002 0.192 -0.435 -0.001 0.429 -0.002 0
UTM_y -0.001 0.146 -0.331 -0.001 0.327 -0.001 0
Random effects:
Name Model
spatial.field SPDE2 model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for the Gaussian observations 352.14 3411.36 5.98e-01 31.27 2357.16 6.06e-01
Range for spatial.field 41377.67 9932.46 2.51e+04 40270.82 64005.40 3.82e+04
Stdev for spatial.field 20.12 3.44 1.43e+01 19.80 27.78 1.91e+01
Deviance Information Criterion (DIC) ...............: 30.47
Deviance Information Criterion (DIC, saturated) ....: 199.46
Effective number of parameters .....................: 90.46
Watanabe-Akaike information criterion (WAIC) ...: 148.63
Effective number of parameters .................: 184.89
Marginal log-Likelihood: 1471.01
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)')plot(output)
#hasil parameter: https://sites.stat.washington.edu/peter/591/Lindgren.pdf
results = inla.spde2.result(output, "spatial.field", Earthquakes_spde)
plot(results[["marginals.range.nominal"]][[1]], type = "l",main = "Nominal range, posterior density")
plot(results[["marginals.variance.nominal"]][[1]], type = "l",main = "Nominal variance, posterior density")
#observe the plots for hyper parameters
list_marginals <- list(
"b0" = output$marginals.fixed$Intercept,
"precision Gaussian obs" = output$marginals.hyperpar$"Precision for the Gaussian observations",
"range" = output$marginals.hyperpar$"Range for spatial.field",
"stdev" = output$marginals.hyperpar$"Stdev for spatial.field")
marginals <- data.frame(do.call(rbind, list_marginals))
marginals$parameter <- rep(names(list_marginals),
times = sapply(list_marginals, nrow))
ggplot(marginals, aes(x = x, y = y)) + geom_line() +
facet_wrap(~parameter, scales = "free") +
labs(x = "", y = "Density") + theme_bw()
Hasil priornya Gaussian obsnya kurang bagus karena tidak berbentuk lonceng, maka di coba2 lagi dengan menetapkan beberapa alternatif untuk gaussian
3.2 Spatiotemporal: interaksinya AR1
formula_ar <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model = Earthquakes_spde, group = spatial.field.group,
control.group = list(model = "ar1"))
output_ar <- inla(formula_ar,data=inla.stack.data(stack, spde=Earthquakes_spde),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack), compute=TRUE))
summary(output_ar)Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data,
quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials =
Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
verbose, ", " lincomb = lincomb, selection = selection, control.compute = control.compute, ", "
control.predictor = control.predictor, control.family = control.family, ", " control.inla =
control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert =
control.expert, ", " control.hazard = control.hazard, control.lincomb = control.lincomb, ", "
control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso =
control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg =
inla.arg, num.threads = num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug = debug, .parent.frame =
.parent.frame)" )
Time used:
Pre = 1.06, Running = 20, Post = 6.15, Total = 27.2
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
Intercept 1.511 0.013 1.486 1.511 1.537 1.511 0
UTM_x 0.000 0.000 0.000 0.000 0.000 0.000 0
UTM_y 0.000 0.000 0.000 0.000 0.000 0.000 0
Random effects:
Name Model
spatial.field SPDE2 model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for the Gaussian observations 156.729 21.559 117.993 155.466 202.699 153.345
Range for spatial.field 576.475 538.334 103.134 420.966 2006.429 240.638
Stdev for spatial.field 0.281 0.334 0.019 0.177 1.160 0.050
GroupRho for spatial.field 0.001 0.701 -0.990 0.006 0.986 -0.999
Deviance Information Criterion (DIC) ...............: -232.34
Deviance Information Criterion (DIC, saturated) ....: 117.89
Effective number of parameters .....................: 6.42
Watanabe-Akaike information criterion (WAIC) ...: -232.88
Effective number of parameters .................: 5.51
Marginal log-Likelihood: 76.97
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)')Nilai DIC RW1: 30.47 dan AR1: -232.34. Nilai DIC yang bagus adalah yang tidak negatif
plot(output_ar)
#hasil parameter: https://sites.stat.washington.edu/peter/591/Lindgren.pdf
result_ar = inla.spde2.result(output_ar, "spatial.field", Earthquakes_spde)
plot(result_ar[["marginals.range.nominal"]][[1]], type = "l",main = "Nominal range, posterior density")
plot(result_ar[["marginals.variance.nominal"]][[1]], type = "l",main = "Nominal variance, posterior density")
list_marginals_ar <- list(
"b0" = output_ar$marginals.fixed$Intercept,
"precision Gaussian obs" = output_ar$marginals.hyperpar$"Precision for the Gaussian observations",
"range" = output_ar$marginals.hyperpar$"Range for spatial.field",
"stdev" = output_ar$marginals.hyperpar$"Stdev for spatial.field",
"rho" = output_ar$marginals.hyperpar$"GroupRho for spatial.field")
marginals_ar <- data.frame(do.call(rbind, list_marginals_ar))
marginals_ar$parameter <- rep(names(list_marginals_ar),
times = sapply(list_marginals_ar, nrow))
ggplot(marginals_ar, aes(x = x, y = y)) + geom_line() +
facet_wrap(~parameter, scales = "free") +
labs(x = "", y = "Density") + theme_bw()
Hasil Pemetaan
#================ Earthquake Prediction ==============================
index_pred <- inla.stack.index(stack,"pred")$data
index_est <- inla.stack.index(stack,"est")$data
Earthquakes_Pred<-data.frame(x=All.Grid1[,1],y=All.Grid1[,2],Time=GroupTime1)
Earthquakes_Est<-data.frame(x=Earthquakes20$UTM_x,y=Earthquakes20$UTM_y, Time=Earthquakes20$Month)
Earthquakes_Pred$pred_mean <- exp(output_ar$summary.fitted.values[index_pred, "mean"])
Earthquakes_Pred$pred_ll <- exp(output_ar$summary.fitted.values[index_pred, "0.025quant"])
Earthquakes_Pred$pred_ul <- exp(output_ar$summary.fitted.values[index_pred, "0.975quant"])
summary(Earthquakes_Pred$pred_mean)
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 4.182 4.462 4.679 4.695 4.918 5.246
Earthquakes_Est$pred_mean <- exp(output_ar$summary.fitted.values[index_est, "mean"])
Earthquakes_Est$pred_ll <- exp(output_ar$summary.fitted.values[index_est, "0.025quant"])
Earthquakes_Est$pred_ul <- exp(output_ar$summary.fitted.values[index_est, "0.975quant"])
summary(Earthquakes_Est$pred_mean)
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 4.332 4.634 4.783 4.761 4.899 5.111
Earthquakes_Values<-rbind(Earthquakes_Est,Earthquakes_Pred)
#head(Earthquakes_Pred)
write.csv(Earthquakes_Pred,"Earthquakes_Pred_AR1.csv")
write_xlsx(Earthquakes_Pred,"Earthquakes_Pred_AR1.xlsx")
Earthquakes_Pred$title <- "By Month"
ggplot() +
geom_tile(data = Earthquakes_Pred, aes(x = x, y = y, fill = pred_mean)) +
scale_fill_gradientn(colours = c("gray99", "yellow", "red")) +
facet_wrap(~ Time, labeller = as_labeller(Month_names), ncol = 4) +
theme_bw() +
ylab("") +
xlab("") +
geom_sf(data = Sumatra_sf_utm1, fill = NA, size=0.01) +
theme(legend.position = "bottom", text = element_text(size = 19)) +
guides(fill = guide_colorbar(title.position = "left", title.vjust = 1,
frame.colour = "black",
barwidth = 20,
barheight = 1.5)) +
theme(panel.background = element_rect(fill = "white", color = NA),
panel.grid.major = element_line(color = "grey95", linetype = "solid"),
panel.grid.minor = element_line(color = "grey95", linetype = "solid"),
axis.text.x = element_blank(), #remove x axis labels
axis.ticks.x = element_blank(), #remove x axis ticks
axis.text.y = element_blank(), #remove y axis labels
axis.ticks.y = element_blank() #remove y axis ticks
) +
labs(fill = "Earthquakes \nMagnitude")
Warna yang terlihat pada peta hasil interpolasi ini memang bisa memberikan kesan oversmoothing, terutama jika perubahan magnitudo gempa bumi terlihat terlalu halus dan tidak mencerminkan pola spasial yang tajam atau abrupt yang mungkin sebenarnya ada di data.
3.3 Spatiotemporal: time: RW1 dan spatial: AR1
hyper.prec <- list(theta = list(prior = "pc.prec", param = c(1, 0.1)))
formula_st <- MAG ~ -1 + Intercept +
f(Time, model="rw1", scale.model = TRUE, hyper =hyper.prec,cyclic=TRUE)+
f(spatial.field, model = Earthquakes_spde, group = spatial.field.group, control.group = list(model = "ar1"))
output_st <- inla(formula_st,data=inla.stack.data(stack, spde=Earthquakes_spde),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack), compute=TRUE))
summary(output_st)Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, quantiles =
quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials = Ntrials, strata =
strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = verbose, ", " lincomb =
lincomb, selection = selection, control.compute = control.compute, ", " control.predictor =
control.predictor, control.family = control.family, ", " control.inla = control.inla, control.fixed =
control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " control.hazard =
control.hazard, control.lincomb = control.lincomb, ", " control.update = control.update, control.lp.scale =
control.lp.scale, ", " control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call
= inla.call, inla.arg = inla.arg, num.threads = num.threads, ", " keep = keep, working.directory =
working.directory, silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug = debug, .parent.frame
= .parent.frame)" )
Time used:
Pre = 1.11, Running = 26.6, Post = 7.04, Total = 34.8
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
Intercept 1.557 0.009 1.539 1.557 1.574 1.557 0
Random effects:
Name Model
Time RW1 model
spatial.field SPDE2 model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for the Gaussian observations 132.54 1.86e+01 99.764 131.264 1.73e+02 128.751
Precision for Time 31060.62 1.38e+05 466.547 7040.024 2.08e+05 1032.240
Range for spatial.field 633.39 6.26e+02 100.642 450.660 2.29e+03 243.136
Stdev for spatial.field 0.39 4.64e-01 0.029 0.247 1.61e+00 0.078
GroupRho for spatial.field 0.00 7.04e-01 -0.989 0.000 9.89e-01 -0.998
Deviance Information Criterion (DIC) ...............: -213.16
Deviance Information Criterion (DIC, saturated) ....: 118.22
Effective number of parameters .....................: 7.31
Watanabe-Akaike information criterion (WAIC) ...: -213.23
Effective number of parameters .................: 6.69
Marginal log-Likelihood: 87.94
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)')plot(output_st)
#hasil parameter: https://sites.stat.washington.edu/peter/591/Lindgren.pdf
result_st = inla.spde2.result(output_st, "spatial.field", Earthquakes_spde)
plot(result_st[["marginals.range.nominal"]][[1]], type = "l",main = "Nominal range, posterior density")
plot(result_st[["marginals.variance.nominal"]][[1]], type = "l",main = "Nominal variance, posterior density")
list_marginals_st <- list(
"b0" = output_st$marginals.fixed$Intercept,
"precision Gaussian obs" = output_st$marginals.hyperpar$"Precision for the Gaussian observations",
"precision for Time" = output_st$marginals.hyperpar$"Precision for Time",
"range" = output_st$marginals.hyperpar$"Range for spatial.field",
"stdev" = output_st$marginals.hyperpar$"Stdev for spatial.field",
"rho" = output_st$marginals.hyperpar$"GroupRho for spatial.field")
marginals_st <- data.frame(do.call(rbind, list_marginals_st))
marginals_st$parameter <- rep(names(list_marginals_st),
times = sapply(list_marginals_st, nrow))
ggplot(marginals_st, aes(x = x, y = y)) + geom_line() +
facet_wrap(~parameter, scales = "free") +
labs(x = "", y = "Density") + theme_bw()
Earthquakes_Pred_st<-data.frame(x=All.Grid1[,1],y=All.Grid1[,2],Time=GroupTime1)
Earthquakes_Est_st<-data.frame(x=Earthquakes20$UTM_x,y=Earthquakes20$UTM_y, Time=Earthquakes20$Month)
Earthquakes_Pred_st$pred_mean <- exp(output_st$summary.fitted.values[index_pred, "mean"])
Earthquakes_Pred_st$pred_ll <- exp(output_st$summary.fitted.values[index_pred, "0.025quant"])
Earthquakes_Pred_st$pred_ul <- exp(output_st$summary.fitted.values[index_pred, "0.975quant"])
summary(Earthquakes_Pred_st$pred_mean)
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 4.702 4.724 4.740 4.744 4.756 4.856
Earthquakes_Est_st$pred_mean <- exp(output_st$summary.fitted.values[index_est, "mean"])
Earthquakes_Est_st$pred_ll <- exp(output_st$summary.fitted.values[index_est, "0.025quant"])
Earthquakes_Est_st$pred_ul <- exp(output_st$summary.fitted.values[index_est, "0.975quant"])
summary(Earthquakes_Est_st$pred_mean)
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 4.700 4.724 4.740 4.746 4.756 5.076
write.csv(Earthquakes_Pred_st,"Earthquakes_Pred_ST.csv")
write_xlsx(Earthquakes_Pred_st,"Earthquakes_Pred_ST.xlsx")
Earthquakes_Pred_st$title <- "By Month"
ggplot() +
geom_tile(data = Earthquakes_Pred_st, aes(x = x, y = y, fill = pred_mean)) +
scale_fill_gradientn(colours = c("gray99", "yellow", "red")) +
facet_wrap(~ Time, labeller = as_labeller(Month_names), ncol = 4) +
theme_bw() +
ylab("") +
xlab("") +
geom_sf(data = Sumatra_sf_utm1, fill = NA, size=0.01) +
theme(legend.position = "bottom", text = element_text(size = 19)) +
guides(fill = guide_colorbar(title.position = "left", title.vjust = 1,
frame.colour = "black",
barwidth = 20,
barheight = 1.5)) +
theme(panel.background = element_rect(fill = "white", color = NA),
panel.grid.major = element_line(color = "grey95", linetype = "solid"),
panel.grid.minor = element_line(color = "grey95", linetype = "solid"),
axis.text.x = element_blank(), #remove x axis labels
axis.ticks.x = element_blank(), #remove x axis ticks
axis.text.y = element_blank(), #remove y axis labels
axis.ticks.y = element_blank() #remove y axis ticks
) +
labs(fill = "Earthquakes \nMagnitude")
3.4 Spatial
#===== 4. spasial: https://punama.github.io/BDI_INLA/
formulas<- MAG ~ -1 + Intercept + UTM_x + UTM_y + f(spatial.field, model=Earthquakes_spde)
ouputs<-inla(formulas, #the formula
data=inla.stack.data(stack,spde=Earthquakes_spde), #the data stack
family= 'gaussian', #which family the data comes from
Ntrials = n, #this is specific to binomial as we need to tell it the number of examined
control.predictor=list(A=inla.stack.A(stack),compute=TRUE), #compute gives you the marginals of the linear predictor
control.compute = list(dic = TRUE, waic = TRUE, config = TRUE), #model diagnostics and config = TRUE gives you the GMRF
verbose = FALSE)
summary(ouputs)Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data,
quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials =
Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
verbose, ", " lincomb = lincomb, selection = selection, control.compute = control.compute, ", "
control.predictor = control.predictor, control.family = control.family, ", " control.inla =
control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert =
control.expert, ", " control.hazard = control.hazard, control.lincomb = control.lincomb, ", "
control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso =
control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg =
inla.arg, num.threads = num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug = debug, .parent.frame =
.parent.frame)" )
Time used:
Pre = 0.99, Running = 5.68, Post = 0.542, Total = 7.22
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
Intercept 1.511 0.013 1.486 1.511 1.537 1.511 0
UTM_x 0.000 0.000 0.000 0.000 0.000 0.000 0
UTM_y 0.000 0.000 0.000 0.000 0.000 0.000 0
Random effects:
Name Model
spatial.field SPDE2 model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for the Gaussian observations 156.722 21.555 118.030 155.448 202.72 153.283
Range for spatial.field 582.742 548.040 99.470 424.221 2037.76 237.808
Stdev for spatial.field 0.285 0.355 0.019 0.175 1.21 0.051
Deviance Information Criterion (DIC) ...............: -232.33
Deviance Information Criterion (DIC, saturated) ....: 117.89
Effective number of parameters .....................: 6.45
Watanabe-Akaike information criterion (WAIC) ...: -232.85
Effective number of parameters .................: 5.54
Marginal log-Likelihood: 77.00
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)')list_marginalss <- list(
"b0" = ouputs$marginals.fixed$Intercept,
"precision Gaussian obs" = ouputs$marginals.hyperpar$"Precision for the Gaussian observations",
"range" = ouputs$marginals.hyperpar$"Range for spatial.field",
"stdev" = ouputs$marginals.hyperpar$"Stdev for spatial.field")
marginalss <- data.frame(do.call(rbind, list_marginalss))
marginalss$parameter <- rep(names(list_marginalss),
times = sapply(list_marginalss, nrow))
ggplot(marginalss, aes(x = x, y = y)) + geom_line() +
facet_wrap(~parameter, scales = "free") +
labs(x = "", y = "Density") + theme_bw()
index_pred <- inla.stack.index(stack,"pred")$data
index_est <- inla.stack.index(stack,"est")$data
Earthquakes_Preds<-data.frame(x=All.Grid1[,1],y=All.Grid1[,2],Time=GroupTime1)
Earthquakes_Est<-data.frame(x=Earthquakes20$UTM_x,y=Earthquakes20$UTM_y, Time=Earthquakes20$Month)
Earthquakes_Preds$pred_mean <- exp(ouputs$summary.fitted.values[index_pred, "mean"])
Earthquakes_Preds$pred_ll <- exp(ouputs$summary.fitted.values[index_pred, "0.025quant"])
Earthquakes_Preds$pred_ul <- exp(ouputs$summary.fitted.values[index_pred, "0.975quant"])
summary(Earthquakes_Preds$pred_mean)
ggplot() +
geom_tile(data = Earthquakes_Preds, aes(x = x, y = y, fill = pred_mean)) +
scale_fill_gradientn(colours = c("gray99", "yellow", "red")) +
theme_bw() +
ylab("") +
xlab("") +
geom_sf(data = Sumatra_sf_utm1, fill = NA, size=0.01) +
theme(legend.position = "bottom", text = element_text(size = 19)) +
guides(fill = guide_colorbar(title.position = "left", title.vjust = 1,
frame.colour = "black",
barwidth = 20,
barheight = 1.5)) +
theme(panel.background = element_rect(fill = "white", color = NA),
panel.grid.major = element_line(color = "grey95", linetype = "solid"),
panel.grid.minor = element_line(color = "grey95", linetype = "solid"),
axis.text.x = element_blank(), #remove x axis labels
axis.ticks.x = element_blank(), #remove x axis ticks
axis.text.y = element_blank(), #remove y axis labels
axis.ticks.y = element_blank() #remove y axis ticks
) +
labs(fill = "Earthquakes \nMagnitude")
3.5 Spatial data
Earthquakes<-read.csv("INLA - Gempa.csv",sep=";")
head(Earthquakes)
#======= Data 2020
Earthquakes20=Earthquakes[Earthquakes$Tahun==2020,]
head(Earthquakes20)
str(Earthquakes20)
#======= Data Preprocessing
Earthquakes20$LAT <- gsub(",", ".", Earthquakes20$LAT)
Earthquakes20$LAT <- as.numeric(Earthquakes20$LAT)
Earthquakes20$LON <- gsub(",", ".", Earthquakes20$LON)
Earthquakes20$LON <- as.numeric(Earthquakes20$LON)
Earthquakes20$MAG <- gsub(",", ".", Earthquakes20$MAG)
Earthquakes20$MAG <- as.numeric(Earthquakes20$MAG)
colnames(Earthquakes20)[colnames(Earthquakes20) == "Bulan"] <- "Month"
str(Earthquakes20)
head(Earthquakes20)
Observed<-data.frame(X=Earthquakes20$LON,Y=Earthquakes20$LAT)
coordinates(Observed) <- ~X+Y # Asumsikan kolom X dan Y adalah kolom lon dan lat
proj4string(Observed) <- CRS("+proj=longlat +datum=WGS84") # Set CRS ke WGS84 (EPSG: 4326)
Observed_utm <- spTransform(Observed, CRS("+proj=utm +zone=48 +datum=WGS84 +units=m +no_defs"))
Observed_utm1<-coordinates(Observed_utm)
colnames(Observed_utm1)<-c("X","Y")
Earthquakes20$UTM_x<-Observed_utm1[,1]
Earthquakes20$UTM_y<-Observed_utm1[,2]
Earthquakes20$IT<-Earthquakes20$Month
head(Earthquakes20)
mesh3 <- inla.mesh.2d(
loc = cbind(Earthquakes20$UTM_x, Earthquakes20$UTM_y), # UTM coordinates of earthquakes
offset = c(50000, 200000), # Min and max distance to mesh boundary
max.edge = c(100000, 1000000) # Maximum triangle edge length
)
coordinates.allyear<-as.matrix(data.frame(x=Earthquakes20$UTM_x,y=Earthquakes20$UTM_y))
control <- list(
results = list(return.marginals.random = TRUE, return.marginals.predictor=TRUE),
compute = list(hyperpar=TRUE, return.marginals.predictor=TRUE, return.marginals=TRUE, dic=TRUE, mlik = TRUE, cpo = TRUE,
po = TRUE, waic=TRUE))
A<-inla.spde.make.A(mesh=mesh3,loc=as.matrix(coordinates.allyear));dim(A)
spde <- inla.spde2.matern(mesh3, alpha=2)
iset <- inla.spde.make.index(name = "spatial.field", spde$n.spde)
stk <- inla.stack(data=list(MAG = log(Earthquakes20$MAG)), #the response
A=list(A,1), #the A matrix; the 1 is included to make the list(covariates)
effects=list(c(list(Intercept=1), #the Intercept
iset), #the spatial indeX
list(UTM_x = Earthquakes20$UTM_x / 1000,
UTM_y = Earthquakes20$UTM_y / 1000)
), tag='dat')
formulas<- MAG ~ -1 + Intercept + UTM_x + UTM_y + f(spatial.field, model=spde)
output_s<-inla(formulas, #the formula
data=inla.stack.data(stk,spde=spde), #the data stack
family= 'gaussian', #which family the data comes from
Ntrials = n, #this is specific to binomial as we need to tell it the number of examined
control.predictor=list(A=inla.stack.A(stk),compute=TRUE), #compute gives you the marginals of the linear predictor
control.compute = list(dic = TRUE, waic = TRUE, config = TRUE), #model diagnostics and config = TRUE gives you the GMRF
verbose = FALSE)
summary(output_s)Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, quantiles =
quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials = Ntrials, strata =
strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = verbose, ", " lincomb =
lincomb, selection = selection, control.compute = control.compute, ", " control.predictor =
control.predictor, control.family = control.family, ", " control.inla = control.inla, control.fixed =
control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " control.hazard =
control.hazard, control.lincomb = control.lincomb, ", " control.update = control.update,
control.lp.scale = control.lp.scale, ", " control.pardiso = control.pardiso, only.hyperparam =
only.hyperparam, ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = num.threads, ", " keep
= keep, working.directory = working.directory, silent = silent, ", " inla.mode = inla.mode, safe =
FALSE, debug = debug, .parent.frame = .parent.frame)" )
Time used:
Pre = 0.873, Running = 0.887, Post = 0.0427, Total = 1.8
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
Intercept 1.508 0.013 1.483 1.508 1.533 1.508 0
UTM_x 0.000 0.000 0.000 0.000 0.000 0.000 0
UTM_y 0.000 0.000 0.000 0.000 0.000 0.000 0
Random effects:
Name Model
spatial.field SPDE2 model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for the Gaussian observations 14487.54 2.81e+04 319.81 6242.64 79505.10 636.88
Theta1 for spatial.field 10.32 2.01e-01 9.92 10.32 10.71 10.32
Theta2 for spatial.field -9.40 1.35e-01 -9.66 -9.40 -9.13 -9.40
Deviance Information Criterion (DIC) ...............: -519.32
Deviance Information Criterion (DIC, saturated) ....: 231.70
Effective number of parameters .....................: 121.56
Watanabe-Akaike information criterion (WAIC) ...: -522.41
Effective number of parameters .................: 96.58
Marginal log-Likelihood: 82.98
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)')#hasil parameter: https://sites.stat.washington.edu/peter/591/Lindgren.pdf
result_s = inla.spde2.result(output_s, "spatial.field", spde)
plot(result_s[["marginals.range.nominal"]][[1]], type = "l",main = "Nominal range, posterior density")
plot(result_s[["marginals.variance.nominal"]][[1]], type = "l",main = "Nominal variance, posterior density")
list_marginalss <- list(
"b0" = output_s$marginals.fixed$Intercept,
"precision Gaussian obs" = output_s$marginals.hyperpar$"Precision for the Gaussian observations",
"theta1" = output_s$marginals.hyperpar$"Theta1 for spatial.field",
"theta2" = output_s$marginals.hyperpar$"Theta2 for spatial.field")
marginalss <- data.frame(do.call(rbind, list_marginalss))
marginalss$parameter <- rep(names(list_marginalss),
times = sapply(list_marginalss, nrow))
ggplot(marginalss, aes(x = x, y = y)) + geom_line() +
facet_wrap(~parameter, scales = "free") +
labs(x = "", y = "Density") + theme_bw()
4. Precision Gaussian
Grid.CoordUTM1<-as.data.frame(coordsGrid.UTM1)
m1<-nrow(Grid.CoordUTM1)
Grid.CoordUTM1<-as.matrix(Grid.CoordUTM1)
GroupTime1<-rep(c(1:12),each=m1)
All.Grid1<-kronecker(matrix(1, 12, 1), Grid.CoordUTM1)
#============ 7.1.1 PENALIZED COMPLEXITY : 10km, 0.1
prior.median.sd1001 = 0.1; prior.median.range10 = 10000
Earthquakes_spde1001 = inla.spde2.pcmatern(mesh=Earthquakes_mesh,
prior.range = c(prior.median.range10, .5),
prior.sigma = c(prior.median.sd1001, .5))
s_index1001 <- inla.spde.make.index(name="spatial.field",
n.spde=Earthquakes_spde1001$n.spde,
n.group=12)
stack_est1001 <- inla.stack(data=list(MAG=log(Earthquakes20$MAG)),
A=list(1,A_est),
effects=list(data.frame(Intercept = rep(1, nrow(Earthquakes20)),
Time=Earthquakes20$Month, UTM_x=Earthquakes20$UTM_x/1000,
UTM_y=Earthquakes20$UTM_y/1000, DEPTH=Earthquakes20$DEPTH),
s_index1001) , tag="est")
A_pred1 <- inla.spde.make.A(mesh=Earthquakes_mesh,
loc=as.matrix(All.Grid1),
group=GroupTime1, #selected day for prediction
n.group=12)
Cov1<-data.frame(UTM_x=All.Grid1[,1]/1000,UTM_y=All.Grid1[,2]/1000,DEPTH=NA,Time=GroupTime1)
stack_pred1001 <- inla.stack(data=list(MAG=NA),
A=list(1, A_pred1),
effects=list(data.frame(Intercept = rep(1, nrow(All.Grid1)),
Time=GroupTime1, UTM_x=Cov1$UTM_x,
UTM_y=Cov1$UTM_y, DEPTH=Cov1$DEPTH), s_index1001),
tag="pred")
stack1001 <- inla.stack(stack_est1001, stack_pred1001)
formula1001 <- MAG ~ -1 + Intercept + UTM_x + UTM_y +
f(spatial.field, model=Earthquakes_spde1001,group=spatial.field.group,
control.group=list(model="rw1"))4.1 Ga ditentuin priornya
output1001 <- inla(formula1001,
data=inla.stack.data(stack1001, spde=Earthquakes_spde1001),
family="gaussian",control.compute = control$compute,
control.predictor=list(A=inla.stack.A(stack1001), compute=TRUE))
summary(output1001)Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data,
quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials =
Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
verbose, ", " lincomb = lincomb, selection = selection, control.compute = control.compute, ", "
control.predictor = control.predictor, control.family = control.family, ", " control.inla =
control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert =
control.expert, ", " control.hazard = control.hazard, control.lincomb = control.lincomb, ", "
control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso =
control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg =
inla.arg, num.threads = num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug = debug, .parent.frame =
.parent.frame)" )
Time used:
Pre = 1.03, Running = 70.3, Post = 3.95, Total = 75.3
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
Intercept 0.447 21.438 -44.814 0.561 45.399 0.511 0
UTM_x -0.002 0.156 -0.365 -0.001 0.358 -0.002 0
UTM_y -0.001 0.109 -0.251 -0.001 0.245 -0.001 0
Random effects:
Name Model
spatial.field SPDE2 model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for the Gaussian observations 523.63 4964.427 4.23e-01 45.45 3586.37 1.49e-01
Range for spatial.field 52442.16 6825.158 4.03e+04 52001.14 67125.41 5.11e+04
Stdev for spatial.field 7.80 0.688 6.55e+00 7.77 9.26 7.69e+00
Deviance Information Criterion (DIC) ...............: 7.83
Deviance Information Criterion (DIC, saturated) ....: 230.95
Effective number of parameters .....................: 121.95
Watanabe-Akaike information criterion (WAIC) ...: 108.05
Effective number of parameters .................: 191.58
Marginal log-Likelihood: 1406.62
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')plot(output1001)
list_marginals <- list(
"b0" = output1001$marginals.fixed$Intercept,
"precision Gaussian obs" = output1001$marginals.hyperpar$"Precision for the Gaussian observations",
"range" = output1001$marginals.hyperpar$"Range for spatial.field",
"stdev" = output1001$marginals.hyperpar$"Stdev for spatial.field"
)
marginals <- data.frame(do.call(rbind, list_marginals))
marginals$parameter <- rep(names(list_marginals),
times = sapply(list_marginals, nrow))
ggplot(marginals, aes(x = x, y = y)) + geom_line() +
facet_wrap(~parameter, scales = "free") +
labs(x = "", y = "Density") + theme_bw()
4.2 loggamma
output1001_lgamma <- inla(
formula1001,
data = inla.stack.data(stack1001,spde=Earthquakes_spde1001),
family = "gaussian", control.compute = control$compute,
control.family = list(
hyper = list(prec = list(prior = "loggamma", param = c(1, 0.01))) # Prior pada precision Gaussian
), control.predictor = list(A = inla.stack.A(stack1001), compute = TRUE)
)
plot(output1001_lgamma)
list_marginals_lgamma <- list(
"b0" = output1001_lgamma$marginals.fixed$Intercept,
"precision Gaussian obs" = output1001_lgamma$marginals.hyperpar$"Precision for the Gaussian observations",
"range" = output1001_lgamma$marginals.hyperpar$"Range for spatial.field",
"stdev" = output1001_lgamma$marginals.hyperpar$"Stdev for spatial.field"
)
marginals_lgamma <- data.frame(do.call(rbind, list_marginals_lgamma))
marginals_lgamma$parameter <- rep(names(list_marginals_lgamma),
times = sapply(list_marginals_lgamma, nrow))
ggplot(marginals_lgamma, aes(x = x, y = y)) + geom_line() +
facet_wrap(~parameter, scales = "free") +
labs(x = "", y = "Density") + theme_bw()
4.3 Gaussian
output1001_lnormal <- inla(
formula1001,
data = inla.stack.data(stack1001,spde=Earthquakes_spde1001),
family = "gaussian", control.compute = control$compute,
control.family = list(
hyper = list(prec = list(prior = "gaussian", param = c(0, 1))) # Prior pada precision Gaussian
), control.predictor = list(A = inla.stack.A(stack1001), compute = TRUE)
)
#summary(output1001_lnormal)
plot(output1001_lnormal)
list_marginals_lnormal <- list(
"b0" = output1001_lnormal$marginals.fixed$Intercept,
"precision Gaussian obs" = output1001_lnormal$marginals.hyperpar$"Precision for the Gaussian observations",
"range" = output1001_lnormal$marginals.hyperpar$"Range for spatial.field",
"stdev" = output1001_lnormal$marginals.hyperpar$"Stdev for spatial.field"
)
marginals_lnormal <- data.frame(do.call(rbind, list_marginals_lnormal))
marginals_lnormal$parameter <- rep(names(list_marginals_lnormal),
times = sapply(list_marginals_lnormal, nrow))
ggplot(marginals_lnormal, aes(x = x, y = y)) + geom_line() +
facet_wrap(~parameter, scales = "free") +
labs(x = "", y = "Density") + theme_bw()
4.4 Prior pcprec (Penalized Complexity Prior for Precision)
output1001_pc <- inla(
formula1001,
data = inla.stack.data(stack1001,spde=Earthquakes_spde1001),
family = "gaussian", control.compute = control$compute,
control.family = list(
hyper = list(prec = list(prior = "pcprec",param = c(1, 0.01) # Uji beberapa kombinasi
))), control.predictor = list(A = inla.stack.A(stack1001), compute = TRUE)
)
plot(output1001_pc)
list_marginals_pc <- list(
"b0" = output1001_pc$marginals.fixed$Intercept,
"precision Gaussian obs" = output1001_pc$marginals.hyperpar$"Precision for the Gaussian observations",
"range" = output1001_pc$marginals.hyperpar$"Range for spatial.field",
"stdev" = output1001_pc$marginals.hyperpar$"Stdev for spatial.field"
)
marginals_pc <- data.frame(do.call(rbind, list_marginals_pc))
marginals_pc$parameter <- rep(names(list_marginals_pc),
times = sapply(list_marginals_pc, nrow))
ggplot(marginals_pc, aes(x = x, y = y)) + geom_line() +
facet_wrap(~parameter, scales = "free") +
labs(x = "", y = "Density") + theme_bw()