Alternative Approaches for Movement Parameterization

The SPoRC package offers multiple approaches for modeling fish movement between regions, each with different complexity and flexibility trade-offs. This vignette demonstrates how to configure movement parameterization using the Setup_Mod_Movement function with various structural options. In general, there are two primary movement model types, which include:

Unstructured Markov: Discrete movement matrices with optional blocking structures across population, age, year, season, and sex dimensions
Continuous-Time Markov Chain (CTMC): CTMC-based movement with diffusion and preference parameters estimated using formula-based approaches

The following demonstrations use the three region sablefish dataset as a basis (three_rg_sable_data). The initial setup establishes the general model dimensions. Note that n_pop = 1 and n_seas = 1 are specified here, but movement can be configured for multi-population and multi-season models in an analogous fashion.

# Load in packages
library(SPoRC)
data("three_rg_sable_data") # load in data

# setup model dimensions
input_list <- Setup_Mod_Dim(years = seq_along(three_rg_sable_data$years),
                            ages = seq_along(three_rg_sable_data$ages),
                            lens = seq(41,99,2),
                            n_regions = three_rg_sable_data$n_regions,
                            n_sexes = three_rg_sable_data$n_sexes,
                            n_fish_fleets = three_rg_sable_data$n_fish_fleets,
                            n_srv_fleets = three_rg_sable_data$n_srv_fleets,
                            n_seas = 1,
                            n_pop = 1,
                            verbose = TRUE
                            )

Like previous vignettes, the Setup_Mod_Dim() function initializes the model structure with 62 years of data (1960-2021), 30 age classes, 3 regions, 2 sexes, 2 fishery fleets (fixed gear and trawl), and multiple survey fleets.

Unstructured Markov

In the simplest case, movement can be parameterized as an unstructured Markov model (move_type = 0), where movement parameters are constant across model partitions and are estimated as discrete transitions between regions. This results in n_regions × (n_regions - 1) parameters estimated.

# Setup movement
input_list <- Setup_Mod_Movement(
  input_list = input_list,
  do_recruits_move = 0, # recruits don't move
  move_type = 0 # unstructured markov
)
length(unique(input_list$map$move_pars)) # number of parameters estimated

Age Blocks

The unstructured Markov model can be extended to incorporate blocking structures, with parameter sharing within blocks to reduce parameter load. In the following, we specify an unstructured Markov model with 2 age blocks, but constant movement across years and sexes. This results in n_regions × (n_regions - 1) × 2 parameters estimated.

# define age blocks
age_blk <- list(c(1:15), c(16:30))
age_blk

# Setup movement
input_list <- Setup_Mod_Movement(
  input_list = input_list,
  do_recruits_move = 0,
  move_type = 0,
  Movement_ageblk_spec = age_blk,
  Movement_yearblk_spec = "constant",
  Movement_sexblk_spec = "constant"
)

length(unique(input_list$map$move_pars)) # number of parameters estimated

Year Blocks

Year blocks can be specified in a similar fashion. Here, 5 year blocks are specified, resulting in n_regions × (n_regions - 1) × 5 parameters estimated.

# define year blocks
yr_blk <- list(c(1:15), c(16:30), c(31:45), c(46:60), c(61:62))
yr_blk

# Setup movement
input_list <- Setup_Mod_Movement(
  input_list = input_list,
  do_recruits_move = 0,
  move_type = 0,
  Movement_ageblk_spec = "constant",
  Movement_yearblk_spec = yr_blk,
  Movement_sexblk_spec = "constant"
)

length(unique(input_list$map$move_pars)) # number of parameters estimated

Sex Blocks

Sex blocks are specified similarly. The following example shows sex-specific movement, resulting in n_regions × (n_regions - 1) × n_sexes parameters estimated.

# define sex blocks
sx_blk <- as.list(1:2)
sx_blk

# Setup movement
input_list <- Setup_Mod_Movement(
  input_list = input_list,
  do_recruits_move = 0,
  move_type = 0,
  Movement_ageblk_spec = "constant",
  Movement_yearblk_spec = "constant",
  Movement_sexblk_spec = sx_blk
)

Blocks Across All Dimensions

Building on the principles described above, movement blocks can also be defined simultaneously across ages, years, and sexes. The example below specifies two age blocks, five year blocks, and sex-specific movement. This configuration results in n_regions × (n_regions - 1) × 2 × 5 × n_sexes movement parameters being estimated. However, this parameterization is likely excessive and may lead to an unstable model solution. Note that this model uses 1 population and 1 season, but population-specific and season-specific movement rates can be specified in an analogous fashion using Movement_popblk_spec and Movement_seasblk_spec.

# define blocks across all dimensions
age_blk <- list(c(1:15), c(16:30))
yr_blk  <- list(c(1:15), c(16:30), c(31:45), c(46:60), c(61:62))
sx_blk  <- as.list(1:2)

# Setup movement
input_list <- Setup_Mod_Movement(
  input_list = input_list,
  do_recruits_move = 0,
  move_type = 0,
  Movement_ageblk_spec = age_blk,
  Movement_yearblk_spec = yr_blk,
  Movement_sexblk_spec = sx_blk
)

Continuous-Time Markov Chain (CTMC)

One potential approach to reduce model parameterization is to represent movement as a continuous-time Markov chain (CTMC) process, in which transitions are governed by a mechanistic framework composed of diffusion (random dispersal) and taxis (directed preference) components. Unlike the unstructured Markov model, CTMC movement processes are defined in continuous time and are converted to annual movement fractions using the matrix exponential, thereby allowing for sequential transitions among regions.

To implement a CTMC-based movement model, an adjacency matrix must first be defined to specify which regions are connected. In the example below, all regions are assumed to be connected, permitting individuals to move among any spatial strata within a given period. An accompanying data frame is then created to define the covariates associated with movement. This data frame must include columns for pop, regions, years, seas, ages, and sexes, plus any additional covariates to be used in the diffusion or preference formulas. The seas column should include all seasons being modeled; projection year covariate values can also be appended to this data frame when n_proj_yrs_devs > 0.

# as.matrix, since Setup_Mod_Movement wants a plain numeric matrix and igraph
# hands back a sparse one
adjacency <- as.matrix(igraph::as_adjacency_matrix(
  igraph::make_graph(
    ~ 1 - 2,
    2 - 3,
    1 - 3
  )
))

# make ctmc data, must include pop, regions, years, seas, ages, sexes columns
ctmc_data <- expand.grid(
  pop     = 1,
  regions = 1:three_rg_sable_data$n_regions,
  years   = seq_along(three_rg_sable_data$years),
  seas    = 1,
  ages    = seq_along(three_rg_sable_data$ages),
  sexes   = 1:three_rg_sable_data$n_sexes
)

Constant Movement

In the code chunk below, movement is defined as arising from a purely diffusive process, which effectively represents constant movement across regions. This formulation results in the estimation of n_regions parameters, providing a more parsimonious alternative to the unstructured Markov approach that estimates n_regions × (n_regions - 1) parameters. Movement among adjacent areas is determined by the defined adjacency matrix, and no directional preference is specified in this case.

The argument area_r = rep(1, 3) specifies how diffusive processes scale with area size. Here, all areas are assumed to be equal. However, when areas differ in size, area_r should be defined as proportional to area, such that smaller areas are associated with higher diffusion rates.

Note that a region-varying diffusion_formula such as ~0 + factor(regions) still redistributes abundance, even with preference_formula = ~0. An edge of the generator holds the diffusion rate of the region a fish is leaving, so the two directions across a boundary differ whenever the regions differ in diffusion, and at equilibrium density settles proportional to 1 / diffusion: abundance accumulates where fish move slowly. That is a real behavioral mechanism, but it does mean the estimated diffusion field and the estimated preference field both act on the equilibrium spatial distribution, so a region-varying diffusion_formula and a region-varying preference_formula are hard to separate from the spatial pattern alone. They differ in how fast that pattern re-establishes, not in the pattern itself.

# constant diffusion, no preference
diffusion_formula  <- ~0 + factor(regions)
preference_formula <- ~0

# Setup movement
input_list <- Setup_Mod_Movement(
  input_list = input_list,
  do_recruits_move = 0,
  move_type = 1,
  ctmc_move_dat = ctmc_data,
  adjacency_mat = adjacency,
  area_r = rep(1, 3),
  diffusion_formula = diffusion_formula,
  preference_formula = preference_formula
)

length(input_list$par$log_move_diffusion_pars)
length(input_list$par$move_preference_pars)

Age-Varying

Movement can also be specified to vary across model partitions. The examples below illustrate how age-varying movement can be represented within the CTMC framework using both linear and spline-based relationships.

Linear

In this example, diffusion is assumed constant across ages, while preference varies linearly by age within each region. This specification results in the estimation of a single diffusion parameter, along with n_regions additional parameters describing age-specific movement preferences.

# constant diffusion, linear age preference
diffusion_formula  <- ~1
preference_formula <- ~0 + factor(regions):ages

# Setup movement
input_list <- Setup_Mod_Movement(
  input_list = input_list,
  do_recruits_move = 0,
  move_type = 1,
  ctmc_move_dat = ctmc_data,
  adjacency_mat = adjacency,
  area_r = rep(1, 3),
  diffusion_formula = diffusion_formula,
  preference_formula = preference_formula
)

length(input_list$par$log_move_diffusion_pars)
length(input_list$par$move_preference_pars)

Spline

In some cases, age-specific movement patterns may be more complex than a simple linear relationship can represent. To capture nonlinear variation, spline-based age-specific movement can be modeled using a spline basis. In the example below, one diffusion parameter is estimated, along with n_regions × 4 parameters describing age-specific movement preferences.

# constant diffusion, spline-based age preference
diffusion_formula  <- ~1
preference_formula <- ~0 + factor(regions):splines2::bSpline(ages, df = 4, intercept = TRUE)

# Setup movement
input_list <- Setup_Mod_Movement(
  input_list = input_list,
  do_recruits_move = 0,
  move_type = 1,
  ctmc_move_dat = ctmc_data,
  adjacency_mat = adjacency,
  area_r = rep(1, 3),
  diffusion_formula = diffusion_formula,
  preference_formula = preference_formula
)

length(input_list$par$log_move_diffusion_pars)
length(input_list$par$move_preference_pars)

Movement Across All Dimensions

Movement can also be specified to vary continuously across all model dimensions. In this example, movement is modeled as a function of region, age, year, and sex. For multi-population models, factor(pop) can be added to the formula in an analogous fashion. This specification results in the estimation of n_regions × 4 × 6 × n_sexes preference parameters.

# spline-based age and year preference, sex-specific
diffusion_formula  <- ~1
preference_formula <- ~0 + factor(regions):
  splines2::bSpline(ages, df = 4, intercept = TRUE):
  splines2::bSpline(years, df = 6, intercept = TRUE):
  factor(sexes)

# Setup movement
input_list <- Setup_Mod_Movement(
  input_list = input_list,
  do_recruits_move = 0,
  move_type = 1,
  ctmc_move_dat = ctmc_data,
  adjacency_mat = adjacency,
  area_r = rep(1, 3),
  diffusion_formula = diffusion_formula,
  preference_formula = preference_formula
)

length(input_list$par$log_move_diffusion_pars)
length(input_list$par$move_preference_pars)

Process Error

Lastly, process error deviations can also be incorporated into movement estimates. In the example below, movement is modeled using a CTMC framework, although process error can similarly be applied to an unstructured Markov model. Movement is specified to vary smoothly across ages using a spline function, while allowing first order autoregressive deviations across years for each region (move_year_re = 'ar1') that are shared across ages (move_age_re = 'none'). Each region’s deviations are one series with its own sd and year correlation in move_pe_pars, shared by both sexes. Three more switches, move_pop_re, move_seas_re and move_sex_re, say how the deviations vary across populations, seasons and sexes, and move_re_years, move_re_ages and their counterparts restrict them to chosen levels. The options are described in the function documentation (?Setup_Mod_Movement); the sharing they allow is shown after the example.

Users may treat the deviations as random effects (integrated out via the Laplace approximation) by naming move_devs in random when calling fit_model(), or leave them out of random as penalized likelihood terms. The sd and correlations are estimated whenever deviations exist; to fix them, map move_pe_pars off in input_list$map and supply the values through input_list$par$move_pe_pars.

Where a deviation sits depends on the movement type. Under an unstructured Markov model it is an offset on the logit of one origin-destination pair, so deviations are estimated for destination regions only (i.e., n_regions - 1), the source region having no deviation term. Under CTMC movement a deviation is added to a region’s own habitat preference, so there is one per region and move_devs holds a single destination level. That one deviation raises the rate into its region from every neighbor and lowers the rate out of it, which is a year in which the region was simply a better place to be, rather than a year in which one boundary happened to be busier. Only differences in preference reach the generator, so adding the same deviation to every region leaves movement unchanged and the level of the field is held down by the process error penalty alone.

Note that because preference enters the generator additively, a deviation larger than an edge’s diffusion rate can potentially drive that rate negative, and with ctmc_diffusion_bounds = "none" the resulting movement fractions end up departing from the [0, 1] while still summing to one (so abundance is not conserved). It is recommended to use a bounded form ("upwind" or "softplus") whenever CTMC deviations are estimated.

# spline age preference with iid year deviations
diffusion_formula  <- ~1
preference_formula <- ~0 + factor(regions):splines2::bSpline(ages, df = 4, intercept = TRUE)

# Setup movement
input_list <- Setup_Mod_Movement(
  input_list = input_list,
  do_recruits_move = 0,
  move_type = 1,
  ctmc_move_dat = ctmc_data,
  adjacency_mat = adjacency,
  area_r = rep(1, 3),
  diffusion_formula = diffusion_formula,
  preference_formula = preference_formula,
  move_year_re = 'ar1',
  move_age_re = 'none'
)

length(input_list$par$log_move_diffusion_pars)
length(input_list$par$move_preference_pars)
length(unique(input_list$map$move_devs))

Every region (origin-destination pair under the unstructured model) has its own deviations. Each of move_pop_re, move_seas_re and move_sex_re takes "none", "iid" (every level its own, independent), "us" (every level its own, with an unstructured correlation across them), or a list of blocks, where the levels in a block share one deviation and the blocks are independent. "none" and "iid" are the two ends of that list, one block of everything and a block per level. Years and ages take no block list, a block on the mean (Movement_yearblk_spec, Movement_ageblk_spec) being where that goes, and the setup refuses such a block when the matching switch is on since the two describe the same thing twice.

move_pe_spec shares the process error across pairs the same way: "est_all" gives every region (pair) its own sd and AR1 correlations, "est_shared" one set for all, and a block list over the pair table one set per block. The pair table is input_list$data$move_pairs, origins and destinations in origin-then-destination order over the edges of the adjacency matrix, which under the CTMC is one row per region.

A block list applies to every region alike, so sharing across sexes cannot differ by region; what can differ by region is the sd. In the three region model below, the first form has the sexes share one series per region, with regions 1 and 2 sharing an sd and year correlation and region 3 holding its own. The second gives every region its own series for each sex, the two sexes correlated within it, with one sd and correlation for all.

# the pair table under the ctmc: one row per region
input_list$data$move_pairs

# the sexes share one series per region; regions 1 and 2 share an sd and correlation, region 3 has its own
input_list <- Setup_Mod_Movement(
  input_list = input_list,
  do_recruits_move = 0,
  move_type = 1,
  ctmc_move_dat = ctmc_data,
  adjacency_mat = adjacency,
  area_r = rep(1, 3),
  diffusion_formula = diffusion_formula,
  preference_formula = preference_formula,
  ctmc_diffusion_bounds = "upwind",
  move_year_re = "ar1",
  move_sex_re = list(1:2), # one block of both sexes, the same as "none"
  move_pe_spec = list(1:2, 3) # regions 1 and 2 share an sd and correlation, region 3 has its own
)
length(levels(input_list$map$move_devs)) # three series of years, one per region
length(levels(input_list$map$move_pe_pars)) # two log sds and two year correlations

# every region its own series, the two sexes correlated within it, one sd and correlation for all
input_list <- Setup_Mod_Movement(
  input_list = input_list,
  do_recruits_move = 0,
  move_type = 1,
  ctmc_move_dat = ctmc_data,
  adjacency_mat = adjacency,
  area_r = rep(1, 3),
  diffusion_formula = diffusion_formula,
  preference_formula = preference_formula,
  ctmc_diffusion_bounds = "upwind",
  move_year_re = "ar1",
  move_sex_re = "us",
  move_pe_spec = "est_shared"
)
length(levels(input_list$map$move_sex_corr_pars)) # one correlation between the sexes, shared by every region