Simulating covariates and outcomes

Jakub Nowosad

2026-03-27

Let’s start by attaching the necessary packages and setting the seed.

library(simsam)
library(terra)
#> terra 1.9.11
set.seed(2024-12-13)

Next, we specify a raster template for the simulations.

rast_grid = terra::rast(xmin = 0, xmax = 200, ymin = 0, ymax = 200, 
                        ncols = 200, nrows = 200)

Simulate

We can simulate any number of covariates using the sim_covariates() function. This function is flexible: it accepts simulation engines created by functions like simulate_gaussian() or simulate_random(). This allows for reusable simulation configurations.

The first example uses simulate_gaussian() with just a range parameter to simulate four covariates for the provided raster grid. The rest of the arguments (e.g., beta representing the average value of the field, nmax representing the number of nearest observations used for kriging simulation, and psill representing the sill of the variogram model) are set to default values.

covariates_one = sim_covariates(rast_grid, n = 4, 
                                method = simulate_gaussian(range = 10))
terra::panel(covariates_one, nr = 1, col.main = "white")

The second example uses a complete variogram model to simulate four covariates for the provided raster grid. Here, we specify not only the range but also the model type (exponential) and the sill of the variogram model.

covariates_two = sim_covariates(rast_grid, n = 4, 
                                method = simulate_gaussian(vgm = gstat::vgm(model = "Exp", psill = 1, range = 10)))
terra::panel(covariates_two, nr = 1, col.main = "white")

The third example uses a nugget effect to simulate four covariates without any spatial correlation.

covariates_three = sim_covariates(rast_grid, n = 4, 
                                  method = simulate_gaussian(vgm = gstat::vgm(model = "Nug", psill = 2, range = 0)))
terra::panel(covariates_three, nr = 1, col.main = "white")

The last example creates a single simulation of a binary outcome (indicators = TRUE) based on a specified variogram model and a beta value. When indicators = TRUE, the simulation produces categorical values (0/1) rather than continuous values. The result depends on the specified variogram model and the beta value (which represents the mean of the underlying Gaussian field).

covariates_four = sim_covariates(rast_grid, n = 1, 
                                method = simulate_gaussian(vgm = gstat::vgm(model = "Exp", psill = 10, range = 100), beta = 30, indicators = TRUE))
terra::plot(covariates_four, nr = 1, col.main = "white")

Blend

The sim_covariates() creates new rasters. Now, we can blend these rasters to create new outcomes with the blend_rasters() function. This function accepts one or more SpatRaster objects, and a formula that specifies how to combine them. Let’s try it out on a few examples.

The first new raster is based on the first set of covariates and a formula that combines them by adding them together.

outcome_one = simsam::blend_rasters(covariates_one, ~ cov1 + cov2 + cov3 + cov4)
terra::plot(outcome_one)

The blend_rasters() function can also be used to further blend the outcomes with additional covariates. Here, we are adding a random noise (cov1 from covariates_three) to the outcome.

outcome_two = simsam::blend_rasters(outcome_one, ~ outcome + cov1, covariates_three)
terra::plot(outcome_two)

The third example shows how to blend the first created outcome with a binary covariate (cov1 from covariates_four). Additionally, it also highlights the flexibility of the formula notation, which allows, for example, the use of various mathematical operations and constants.

outcome_three = simsam::blend_rasters(outcome_one, ~ outcome * (cov1 + 0.2), covariates_four)
terra::plot(outcome_three)

Finally, the last example shows that the blend_rasters() function can be used in a pipe-like manner.

outcome_four = simsam::blend_rasters(covariates_two, ~ cov1 + cov2 * cov3 + cov4) |>
  simsam::blend_rasters(~ outcome * (cov1 + 0.2), covariates_four)
terra::plot(outcome_four)