Install INLA

install.packages("INLA",repos=c(getOption("repos"),INLA="https://inla.r-inla-download.org/R/stable"), dep=TRUE)

Libraries Need

library(sf)
library(sp)
library(INLA)
library(raster)
library(spatstat)
library(dplyr)
library(terra)

Input Data

#This data needs be changed if the study area changes

explosions <- st_read("D:\\Ukraine_GradProject\\Explosion_Model\\Explosion_Model\\Tile_Effects\\exp_model_V4.shp") # Reads in Explosions Point Shapefile

## Reading layer `exp_model_V4' from data source 
##   `D:\Ukraine_GradProject\Explosion_Model\Explosion_Model\Tile_Effects\exp_model_V4.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 84324 features and 44 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 31.3517 ymin: 45.3854 xmax: 39.5732 ymax: 50.6984
## Geodetic CRS:  WGS 84

boundary <- st_read("D:\\Ukraine_GradProject\\Explosion_Model\\Explosion_Model\\Tile_Effects\\explosiontiles.shp") # Reads in Boundary Shapefile

## Reading layer `explosiontiles' from data source 
##   `D:\Ukraine_GradProject\Explosion_Model\Explosion_Model\Tile_Effects\explosiontiles.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 26 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 1586587 ymin: 2026501 xmax: 2096973 ymax: 2661127
## Projected CRS: Europe_Lambert_Conformal_Conic

# Coordinate System in Use: Europe_Lambert_Conformal_Conic (Meters)
explosions <- st_transform(explosions, st_crs(boundary)) # Coverts the CRS into the same CRS as the Boundary file

Data Formation

# Extract Coordinates and Covariates
coords <- st_coordinates(explosions) # Extract the Coordinates from the Explosions
# Create a Dataframe with all the Coordinates and Covariate Information
# Covariate Information needs to have distance to all the covariates, and the population number of the closest city to the explosion point
data <- explosions %>% 
  mutate(x = coords[, 1], y = coords[, 2]) %>%
  select(x, y, DIST_ROAD, DIST_FL, DIST_RAILS, DIST_HOS, DIST_PP, RASTERVALU)

Creating the Spatial Random Effect

# Create Spatial Mesh
mesh <- inla.mesh.2d(loc = coords,
                     max.edge = c(10000,100000),
                     cutoff = 5000)
plot(mesh)

# Stochastic Partial Differential Equation (SPDE) Model
# Creates the Matern Covariance Function for the Spatial Random Effect
spde <- inla.spde2.pcmatern(mesh = mesh, alpha = 2, 
                            prior.range = c(5000, 0.07),  
                            prior.sigma = c(1, 0.1))  
# The Prior Range is larger than 5 km with a 7% Probability it is smaller
# Prior Range includes the distance (5000 m or 5 km) at which spatial correlation drops
# Prior Sigma is that standard deviation is larger than 1 with a probability of 1% (Low Variability in the Spatial Field)

Grid Information for Raster Image and INLA

grid_res <- 250  # Cell Size Resolution of 1000 meters
boundary_extent <- st_bbox(boundary)
x_seq <- seq(boundary_extent["xmin"], boundary_extent["xmax"], by = grid_res)
y_seq <- seq(boundary_extent["ymin"], boundary_extent["ymax"], by = grid_res)
grid_coords <- expand.grid(x = x_seq, y = y_seq)

# Convert to sf object and clip to boundary
# Used in creating final raster image
grid_sf <- st_as_sf(grid_coords, coords = c("x", "y"), crs = st_crs(boundary))
grid_sf <- st_intersection(grid_sf, boundary)

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

# Convert to matrix for INLA
grid_mat <- st_coordinates(grid_sf)

# Drops the Geometry of the Data
data_no_geom <- st_drop_geometry(data)
# Drops the X, Y geometry from the dataframe
covariates_obs <- data_no_geom %>%
  select(-x, -y)

Scaling Covariate Values

# Covariate Values were in Kilometers (Still were too big for INLA to Run). All are divided by the Mean of each covariate besides Population(log+1)
covariates_obs$DIST_ROAD <- (covariates_obs$DIST_ROAD) / mean(covariates_obs$DIST_ROAD)
covariates_obs$DIST_FL <- (covariates_obs$DIST_FL) / mean(covariates_obs$DIST_FL)
covariates_obs$Population <- log(covariates_obs$RASTERVALU + 1)
covariates_obs$DIST_RAILS <- (covariates_obs$DIST_RAILS) / mean(covariates_obs$DIST_RAILS)
covariates_obs$DIST_HOS <- (covariates_obs$DIST_HOS) / mean(covariates_obs$DIST_HOS)
covariates_obs$DIST_PP <- (covariates_obs$DIST_PP) / mean(covariates_obs$DIST_PP)

Projection Matrix for Observed Locations

A_obs <- inla.spde.make.A(mesh = mesh, loc = coords) # Projects the Spatial Random Effect onto the observations

# Index for the spatial random field
spde_index <- inla.spde.make.index(name = "spatial_field", n.spde = spde$n.spde)

# Build observation stack
stack_obs <- inla.stack(
  data = list(y = rep(1, nrow(covariates_obs))),  # Measures Count Data; y=1 for each observation point
  A = list(A_obs, 1), # Projects the Spatial Random Effect onto the observations
  effects = list(
    spatial_field = spde_index$spatial_field, # Links the Spatial Random Effect from the SPDE to the observations
    covariates = covariates_obs # Fixed Effects
  ),
  tag = "obs"
)

Projection Matrix for Prediction Locations

A_pred <- inla.spde.make.A(mesh = mesh, loc = grid_mat)

Prediction Covariates

# This data needs be changed if the study area changes
# These raster layers were created in ArcGIS Pro using the Euclidean Distance tool for the distances and WorldPop data for the Population.
DIST_ROAD <- rast("D:\\Ukraine_GradProject\\Explosion_Model\\Explosion_Model\\Covariates\\roads_ukr_250m.tif")
DIST_FL <- rast("D:\\Ukraine_GradProject\\Explosion_Model\\Explosion_Model\\Covariates\\frontline_ukr_250m.tif")
DIST_RAILS <- rast("D:\\Ukraine_GradProject\\Explosion_Model\\Explosion_Model\\Covariates\\rails_ukr_250m.tif")
DIST_HOS <- rast("D:\\Ukraine_GradProject\\Explosion_Model\\Explosion_Model\\Covariates\\hos_ukr_250m.tif")
DIST_PP <- rast("D:\\Ukraine_GradProject\\Explosion_Model\\Explosion_Model\\Covariates\\pp_ukr_250m.tif")
Population <- rast("D:\\Ukraine_GradProject\\Explosion_Model\\Explosion_Model\\Covariates\\new_pop_v2_250.tif")

# grid_mat is a matrix with columns x and y
coords_pred_vect <- vect(grid_mat, type = "points", crs = crs(DIST_ROAD))

# Extract raster values at prediction points
DIST_ROAD_vals <- raster::extract(DIST_ROAD, coords_pred_vect)[,2]
DIST_FL_vals <- raster::extract(DIST_FL, coords_pred_vect)[,2]
DIST_RAILS_vals <- raster::extract(DIST_RAILS, coords_pred_vect)[,2]
DIST_HOS_vals <- raster::extract(DIST_HOS, coords_pred_vect)[,2]
DIST_PP_vals <- raster::extract(DIST_PP, coords_pred_vect)[,2]
Population_vals <- raster::extract(Population, coords_pred_vect)[,2]

#Covariate Dataframe
covariate_df <- data.frame(
  DIST_ROAD = DIST_ROAD_vals,
  DIST_FL = DIST_FL_vals,
  DIST_RAILS = DIST_RAILS_vals,
  DIST_HOS = DIST_HOS_vals,
  DIST_PP = DIST_PP_vals,
  Population = Population_vals
)

stack_pred <- inla.stack(
  data = list(y = NA),  # NA tells INLA that we are predicting here
  A = list(A_pred, 1), # Projects the Spatial Random Field(on the mesh) to the Prediction Locations
  effects = list(
    spatial_field = spde_index$spatial_field, # Links the Spatial Random Effect from the SPDE to the Predictions
    covariates = covariate_df # Covariate dataframe that has all the prediction values in it
  ),
  tag = "pred"
)

Combine both stacks for the INLA model

stack_full <- inla.stack(stack_obs, stack_pred)

Creating Model Formula and Running INLA

Model Formula

formula <- y ~ DIST_ROAD + DIST_FL + DIST_RAILS + DIST_HOS + DIST_PP + Population + f(spatial_field, model = spde) # No Intercept is Used; All of the Covariates are used and the Spatial Random Effect from SPDE is added

# Model Formula
formula <- y ~ DIST_ROAD + DIST_FL + DIST_RAILS + DIST_HOS + DIST_PP + Population +
  f(spatial_field, model = spde)

# Running INLA
result3 <- inla(
  formula,
  family = "poisson",
  data = inla.stack.data(stack_full),
  control.predictor = list(
    A = inla.stack.A(stack_full),
    compute = TRUE
  ),
  control.inla = list(int.strategy = "eb", h = 0.01),
  control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE),
  verbose = TRUE
)

## Warning in inla.core(formula = formula, family = family, contrasts = contrasts, : The A-matrix in the predictor (see ?control.predictor) is specified
##   but an intercept is in the formula. This will likely result
##   in the intercept being applied multiple times in the model, and is likely
##   not what you want. See ?inla.stack for more information.
##   You can remove the intercept adding ``-1'' to the formula.

Results

result3$summary.fixed # Shows the results of the fixed covariates

##                      mean          sd   0.025quant      0.5quant  0.975quant
## (Intercept)  1.279111e-06 0.006339425 -0.012423766  1.279111e-06 0.012426324
## DIST_ROAD   -9.270921e-06 0.003651255 -0.007165599 -9.270921e-06 0.007147057
## DIST_FL     -6.484392e-06 0.003265063 -0.006405890 -6.484392e-06 0.006392921
## DIST_RAILS  -1.040186e-05 0.003483589 -0.006838112 -1.040186e-05 0.006817308
## DIST_HOS    -1.259424e-06 0.005299950 -0.010388971 -1.259424e-06 0.010386452
## DIST_PP     -1.946742e-05 0.006509315 -0.012777490 -1.946742e-05 0.012738555
## Population  -7.252358e-07 0.001674852 -0.003283375 -7.252358e-07 0.003281924
##                      mode kld
## (Intercept)  1.279111e-06   0
## DIST_ROAD   -9.270921e-06   0
## DIST_FL     -6.484392e-06   0
## DIST_RAILS  -1.040186e-05   0
## DIST_HOS    -1.259424e-06   0
## DIST_PP     -1.946742e-05   0
## Population  -7.252358e-07   0

result3$summary.hyperpar # Shows the results of the priors from the SPDE

##                                 mean           sd   0.025quant     0.5quant
## Range for spatial_field 4.457203e+04 1.137068e+05 2.171664e+03 1.724062e+04
## Stdev for spatial_field 3.376556e-03 2.549372e-03 4.620959e-04 2.706415e-03
##                           0.975quant         mode
## Range for spatial_field 2.622710e+05 5.055867e+03
## Stdev for spatial_field 9.902459e-03 1.358932e-03

## Creating Predictions from the Posterior Distribution (Mean)
index_pred <- inla.stack.index(stack_full, "pred")$data
predictions <- result3$summary.fitted.values[index_pred, "mean"]


# Create Column for Predictions
grid_sf$predictions <- predictions

# Exponentiation of the Explosion Count Intensity
grid_sf$predictions_exp <- exp(grid_sf$predictions)

## Creating the Raster Image
r_pred_exp <- raster(extent(boundary), res = grid_res, crs = st_crs(boundary)$proj4string)
r_pred_exp <- rasterize(grid_sf, r_pred_exp, field = "predictions_exp", fun = mean, na.rm = TRUE)

## Plot the Predicted Explosion Count Intensity Plot
plot(r_pred_exp, main = "Predicted Explosion Intensity Count", col = terrain.colors(100))
plot(st_geometry(boundary), add = TRUE, border = "black", lwd = 2)

ExplosionModel_INLA

Michael Pegg

2026-02-18