Simulating data from the model
Loading required libraries and data
library(JM.landslide) # Contains functions for joint modeling of landslide counts and sizes data.
library(spdep) # Provides functions for spatial dependence analysis.
library(INLA) # Used here for creating the neighborhood structures.
library(ggplot2)
library(gridExtra)
library(viridis)
# The dataset contains the necessary spatial information, such as the shapefile with spatial locations and coordinates.
data("Wenchuan_info_used_for_simulation")Creating the adjacency graph and precision matrix
This below code constructs the precision matrix Q, based
on the spatial dependency structures of the data.
nb2INLA("adjgraph-sim.txt", poly2nb(shp_selected, queen = FALSE, row.names = shp_selected$SU_ID))
adjacensy <- inla.read.graph(filename = "adjgraph-sim.txt")
N <- adjacensy$n
diag.Q <- diag(adjacensy$nnbs, N)
A.Q <- matrix(0, nrow = N, ncol = N)
for (i in 1:N) {
A.Q[i, adjacensy$nbs[[i]]] <- 1
}
Q <- diag.Q - A.Q # Precision matrixHere:
diag.Q: Creates a diagonal matrix where the diagonal elements represent the number of neighbors for each nodeA.Q: Adjacency matrix- The final precision matrix,
Q = diag.Q - A.Q, is obtained by subtracting the adjacency matrix from the diagonal matrix
Setting fixed parameters
Fixed the model hyperparameters \(\kappa_{\boldsymbol w_1}\)
(kappa_w1), \(\kappa_{\boldsymbol
w_2}\) (kappa_w2), \(\kappa_{\boldsymbol \eta}\)
(kappa_eta), \(\kappa_{\boldsymbol \mu }\)
(kappa_mu), the model intercepts \(\gamma_1\) (intercept1) and
\(\gamma_2\) (intercept2),
the sharing correlation parameter \(\beta\) (beta), the regression
coefficients \(\boldsymbol\beta_1\)
(beta1) and \(\boldsymbol\beta_2\) (beta2),
and the parameters of the sub-asymptotic distribution \(\boldsymbol \Theta_A\)
(hyper.mu) to reasonable values. Additionally, simulated
covariate matrices \(\boldsymbol Z_1\)
(Z1) and \(\boldsymbol
Z_2\) (Z2) were generated, which contain the fixed
effects for landslide counts and sizes.
# Set fixed parameters
kappa_w1 <- 10
kappa_w2 <- 5
kappa_eta <- 10
kappa_mu <- 5
intercept1 <- 2
intercept2 <- 4
beta <- 1
other.hyper <- list(
kappa_w1 = kappa_w1,
kappa_w2 = kappa_w2,
kappa_eta = kappa_eta,
kappa_mu = kappa_mu,
intercept1 = intercept1,
intercept2 = intercept2,
beta = beta
)
hyper.mu<- c(20,0.1)
set.seed(1)
Z1 <- mvtnorm::rmvnorm(N, mean = rep(0, 3))
Z2 <- mvtnorm::rmvnorm(N, mean = rep(0, 2))
beta1 <- runif(ncol(Z1), -1, 1)
beta2 <- runif(ncol(Z2), -1, 1)Simulating data
Given a fixed set of parameter values, we can use forward sampling
techniques to simulate data from the model, which is implemented in the
function sim_mark_function. This function simulates
landslide count and size data according to the specified model
introduced in the Introduction section, based on the given precision
matrix, hyperparameters, and covariates. For example, setting
model_type = "jSp" simulates the joint spatial model
described earlier, whereas model_type = "FE" simulates the
data using only the available covariate information, excluding any ICAR
random effects. Additionally, setting mark_dist = "eGPD"
simulates the model where the size distribution \(F_A\) follows an extended generalized
Pareto distribution. Other options include
mark_dist = "bGPD" for a mixture of Beta and generalized
Pareto distributions, and mark_dist = "tgGPD" for a mixture
of truncated Gamma and generalized Pareto distributions. For more
details, see Yadav et
al. (2024).
set.seed(1)
Sim_mark_data <- JM.landslide::sim_mark_function(
Q = Q,
hyper.mu = hyper.mu,
other.hyper = other.hyper,
beta1 = beta1, beta2 = beta2,
Z1 = as.matrix(Z1), Z2 = as.matrix(Z2),
mark_dist = "eGPD",
model_type = "jSp"
)
# Access results
Y <- Sim_mark_data$Y
A <- Sim_mark_data$A
mu <- Sim_mark_data$mu
eta <- Sim_mark_data$eta
W1 <- Sim_mark_data$W1
W2 <- Sim_mark_data$W2Plots
generate_histogram <- function(data, var_name, binwidth = 30, title = NULL, fill_color = "blue", y_label = "Frequency") {
ggplot(data, aes_string(x = var_name)) +
geom_histogram(binwidth = binwidth, fill = fill_color, alpha = 0.7, color = "black") +
labs(title = title, y = y_label) + # Add custom y-axis label
theme_minimal()
}
Sim_data <- data.frame(Y = Sim_mark_data$Y,
A = Sim_mark_data$A,
mu = Sim_mark_data$mu,
eta = Sim_mark_data$eta,
W1 = Sim_mark_data$W1,
W2 = Sim_mark_data$W2)
p1 <- generate_histogram(Sim_data, var_name = "Y", binwidth = 30, title = "Histogram of Y", fill_color = "blue", y_label = "Frequency")
p2 <- generate_histogram(Sim_data, var_name = "A", binwidth = 30, title = "Histogram of A", fill_color = "green", y_label = "Frequency")
p3 <- generate_histogram(Sim_data, var_name = "eta", binwidth = 0.5, title = "Histogram of eta", fill_color = "red", y_label = "Frequency")
p4 <- generate_histogram(Sim_data, var_name = "mu", binwidth = 0.5, title = "Histogram of mu", fill_color = "purple", y_label = "Frequency")
p5 <- generate_histogram(Sim_data, var_name = "W1", binwidth = 0.5, title = "Histogram of W1", fill_color = "orange", y_label = "Frequency")
p6 <- generate_histogram(Sim_data, var_name = "W2", binwidth = 0.5, title = "Histogram of W2", fill_color = "yellow", y_label = "Frequency")
grid.arrange(p1, p2, p3, p4, p5, p6, nrow = 3, ncol = 2)
Here:
Y: Simulated landslide counts.A: Simulated landslide sizes.mu,eta: Latent variables representing the log-linear predictors for counts and sizes, respectively.W1,W2: Spatially structured latent ICAR random effects for counts and sizes, respectively.
