Skip to contents

Specify a first-order random walk temporal random effect. RW1 assumes that the difference between consecutive time points is normally distributed: phi[t] - phi[t-1] ~ N(0, sigma^2).

This is an intrinsic model (improper prior) that captures smooth temporal trends. It is equivalent to an ICAR model on a 1D chain.

Usage

temporal_rw1(
  time_var,
  group_var = NULL,
  cyclic = FALSE,
  shared = TRUE,
  parameterization = c("centered", "noncentered")
)

Arguments

time_var

Name of the time variable in data. Can be a formula (e.g., ~ year) or a character string (e.g., "year").

group_var

Optional name of grouping variable for panel data. If provided, separate random walks are estimated for each group. Each group's walk sums to zero over time, so it carries the group's trend and not its level; add a group intercept such as (1 | site) for the levels.

cyclic

Logical; if TRUE, the random walk wraps around (first and last time points are neighbors). Useful for seasonal data.

shared

Logical; if TRUE (default), temporal effect enters both numerator and denominator linear predictors identically.

parameterization

Coordinate the walk is sampled in. "centered" (default) samples the effects themselves; "noncentered" samples z ~ N(0, I) and reaches the effects through sigma * A^-1 z, where A stacks the walk's difference operator on 1'/sqrt(n), so the variance enters the linear predictor through the transform rather than through a quadratic form. The two describe the same model and differ in the geometry the sampler sees: non-centred helps where the variance is shrunk toward zero and leaves a funnel the centred coordinate cannot descend, and can be worse where the data pin the field. Measured on a 500-observation fixture with no temporal signal, an RW1's bulk ESS for tau_temporal goes 409 to 805, while an RW2 picks up 23 divergences of 2000 for no gain, which is why this is an argument rather than a default. With group_var the prior still augments one global constant, so the transform carries the group contrasts as free directions beside the scaled ones and the group levels stay unshrunk, as they are in the centred coordinate.

Value

A ratiod_temporal object

Details

The RW1 precision matrix Q has the form (for T time points): $$Q[t,t] = 1 \text{ for } t = 1, T$$ $$Q[t,t] = 2 \text{ for } 1 < t < T$$ $$Q[t,t-1] = Q[t-1,t] = -1$$

This is a rank-deficient matrix (rank T-1), so a sum-to-zero constraint is applied for identifiability.

Examples

# Create temporal RW1 specification
temporal_rw1("year")
#> tulpaRatio temporal specification
#> ============================
#> 
#> Type: RW1 (First-order Random Walk) 
#> Time variable: year 
#> Shared: Yes (enters both processes) 
temporal_rw1("month", cyclic = TRUE)
#> tulpaRatio temporal specification
#> ============================
#> 
#> Type: RW1 (First-order Random Walk) 
#> Time variable: month 
#> Cyclic: Yes
#> Shared: Yes (enters both processes) 

# \donttest{
# Simple temporal trend: a random walk shared by count and effort
set.seed(123)
year <- rep(2010:2021, each = 4)
trend <- cumsum(rnorm(12, sd = 0.25))[year - 2009]
x <- rnorm(48)
df <- data.frame(
  year = year,
  x = x,
  count = rpois(48, exp(3 + 0.3 * x + trend)),
  effort = rgamma(48, shape = 5, rate = exp(-trend))
)

fit <- tratio(
  count | effort ~ x,
  data = df,
  family = ratiod_poisson_gamma(),
  temporal = temporal_rw1("year"),
  mode = "hmc",
  control = list(iter = 400, chains = 1, seed = 1)
)

# Panel data: separate trends per site around each site's level
set.seed(124)
site <- rep(paste0("site", 1:5), times = 12)
s <- match(site, paste0("site", 1:5))
t <- rep(1:12, each = 5)
walks <- apply(matrix(rnorm(12 * 5, sd = 0.2), 12, 5), 2, cumsum)
eff <- rnorm(5, sd = 0.5)[s] + walks[cbind(t, s)]
x <- rnorm(60)
df_panel <- data.frame(
  year = 2014 + t,
  site = site,
  x = x,
  count = rpois(60, exp(2.5 + 0.3 * x + eff)),
  effort = rgamma(60, shape = 3, rate = exp(-eff))
)

# The per-site walks carry the trends, (1 | site) the site levels
fit2 <- tratio(
  count | effort ~ x + (1 | site),
  data = df_panel,
  family = ratiod_poisson_gamma(),
  temporal = temporal_rw1("year", group_var = "site"),
  mode = "hmc",
  control = list(iter = 400, chains = 1, seed = 1)
)

# Cyclic for monthly seasonal data
set.seed(125)
month <- rep(1:12, times = 4)
season <- 0.5 * sin(2 * pi * month / 12)
x <- rnorm(48)
df_monthly <- data.frame(
  month = month,
  x = x,
  count = rpois(48, exp(2.5 + 0.3 * x + season)),
  effort = rgamma(48, shape = 2, rate = exp(-season))
)

fit3 <- tratio(
  count | effort ~ x,
  data = df_monthly,
  family = ratiod_poisson_gamma(),
  temporal = temporal_rw1("month", cyclic = TRUE),
  mode = "hmc",
  control = list(iter = 400, chains = 1, seed = 1)
)
# }