fit_dynamic_model: Fit a Bayesian dynamic count / binomial / multinomial...

In DynCount: Bayesian Dynamic Models for Count Time Series

Fit a Bayesian dynamic count / binomial / multinomial time-series model

Description

Fits a GMRF state-space model in which a latent trajectory z_t evolves as either a first-order random walk (latent_dynamics = "rw", the default) or a stationary AR(1) process (latent_dynamics = "ar1"), and the observations are linked to it through a Poisson (log link), binomial (logit link) or multinomial (additive-log-ratio link) observation model. See DynCount-package for an overview.

Usage

fit_dynamic_model(
  y,
  family = c("poisson", "binomial", "multinomial"),
  trials = NULL,
  innovations = c("gaussian", "t", "mixture", "sv"),
  latent_dynamics = c("rw", "ar1"),
  include_mu = FALSE,
  zeros = c("none", "inflated", "missing"),
  zero_inflation = FALSE,
  prior = dynamic_prior(),
  nsave = 4000,
  nburn = 1000,
  thin = 1,
  horizon = 0L,
  forecast_trials = NULL,
  forecast_offset = NULL,
  offset = NULL,
  baseline = "largest",
  verbose = FALSE,
  seed = NULL
)

Arguments

y	For the Poisson and binomial families a numeric vector of non-negative integer observations (counts, or numbers of successes). For the multinomial family an n \times K matrix (or data frame) of non-negative integer category counts with one row per time point; its column names, if any, are used as category labels.
family	Observation model, "poisson" (default), "binomial" or "multinomial".
trials	For family = "binomial", the number of trials. Either a single number (recycled) or a vector the same length as y. Each y must not exceed its number of trials. For the multinomial family the row totals of y are used and this argument must be left NULL; for the Poisson family it is ignored with a warning.
innovations	Distribution of the latent increments, one of "gaussian" (default), "t", "mixture", "sv". See DynCount-package for details. "sv" requires the stochvol package.
latent_dynamics	Latent state evolution: "rw" (default; a first-order random walk, \rho = 1 fixed) or "ar1" (a stationary AR(1) with \rho sampled on (-1, 1)). "ar1" always includes an intercept (include_mu is forced to TRUE).
include_mu	Logical; include a scalar \mu in the state equation z_t = \mu + \rho z_{t-1} + \varepsilon_t. It is a drift under latent_dynamics = "rw" (\rho = 1) and an intercept under "ar1". With FALSE (default) \mu = 0 and is not sampled; with TRUE it has the Gaussian prior in dynamic_prior(). When latent_dynamics = "ar1" this is forced to TRUE regardless of the value supplied: the intercept gives the process a non-zero stationary mean \mu / (1 - \rho), without which the zero-mean stationary assumption is rarely appropriate for a log-rate/logit series.
zeros	Zero handling for the Poisson and binomial families: "none" (default), "inflated" (zero inflation with a time-constant gate-open probability; see structural_zero_prob()), or "missing" (observed zeros treated as missing data). Must be "none" for the multinomial family. See Details.
zero_inflation	A single TRUE or FALSE. TRUE is shorthand for zeros = "inflated"; combining it with a different explicit zeros is an error. Unlike the zero_inflation argument of simulate_dynamic_poisson() and simulate_dynamic_binomial(), it is not a probability.
prior	A dynamic_prior() object giving the prior hyperparameters. For the multinomial family the same prior is applied independently to every non-baseline category.
nsave	Number of posterior draws to keep. Default 4000.
nburn	Number of burn-in iterations. Default 1000.
thin	Thinning interval: one draw is kept every thin iterations, so the sampler runs nburn + nsave * thin iterations in total and retains nsave draws. Default 1.
horizon	Forecast horizon H (a non-negative integer). When H >= 1, an H-step forecast is simulated after sampling and stored for retrieval with forecast.dynamic_fit(); when H = 0 (default) nothing is stored, and forecasts can still be computed later with forecast(fit, horizon = H).
forecast_trials	For the binomial and multinomial families, the number of trials (binomial) or the total count per period (multinomial) over the forecast horizon (length 1, recycled, or length horizon). If omitted it defaults to the last observed number of trials (binomial) or the last non-zero row total (multinomial). Used only when horizon >= 1.
forecast_offset	Known offset over the forecast horizon. For the Poisson and binomial families a vector of length 1 (recycled) or horizon; for the multinomial family a scalar, a length K - 1 vector (one constant per non-baseline category) or an H \times (K - 1) matrix on the ALR scale (columns aligned with the non-baseline categories in their original order). Defaults to 0, with a warning if offset is non-zero. Used only when horizon >= 1.
offset	Optional known per-observation offset on the linear-predictor scale. For the Poisson family (length 1 or n) this is a log-exposure term, so the mean is \exp(\mathrm{offset}_t + z_t); for the binomial family it shifts the logit. For the multinomial family it is a scalar, a length K - 1 vector (one constant per non-baseline category) or an n \times (K - 1) matrix of per-category shifts on the ALR scale (columns aligned with the non-baseline categories in their original order; the baseline has no offset). Default NULL (no offset).
baseline	Multinomial family only: the baseline category, given as a column index or a column name of y, or "largest" (default) to use the category with the largest total count over the series (ties broken by column order). A category with many observations makes the ALR transform numerically well behaved. Note that the model is not invariant to this choice: the latent dynamics are placed on the log-ratios relative to the baseline, so if the baseline's own share moves a lot every ALR series inherits that movement. A large category with a stable share is the natural choice. For the other families a non-default value is ignored with a warning.
verbose	Logical; print a progress bar. Default FALSE.
seed	Optional random seed for reproducibility. The previous state of the global random number generator is restored after fitting.

Details

Estimation. The model is estimated by Metropolis-within-Gibbs MCMC. The latent states are updated one site at a time by adaptive random-walk Metropolis steps that use their Gaussian Markov random field (GMRF) full conditionals. States without an observation (the initial state, structural zeros, zero-total rows) are drawn exactly from their Gaussian full conditionals. The innovation parameters, the drift/intercept \mu and the AR(1) coefficient \rho are updated by Gibbs steps (with a Metropolis step for the Student-t degrees of freedom).

Multinomial family. With family = "multinomial", y is an n \times K matrix of category counts (rows = time, columns = categories, K \ge 2); the row totals N_t are treated as known trials. One category b is the baseline, and the remaining K - 1 categories each get their own latent additive-log-ratio (ALR) series z_{t,k} = \log(p_{t,k} / p_{t,b}), so that

y_t \sim \mathrm{Multinomial}(N_t, p_t), \qquad p_{t,k} = \frac{\exp(o_{t,k} + z_{t,k})}{1 + \sum_{j \ne b} \exp(o_{t,j} + z_{t,j})}, \qquad p_{t,b} = \frac{1}{1 + \sum_{j \ne b} \exp(o_{t,j} + z_{t,j})},

with known offsets o_{t,k} (zero by default). The chosen latent dynamics, innovation structure and drift/intercept setting apply to every ALR series, but no parameters are shared across categories. Each series has its own innovation variance (and, where relevant, degrees of freedom, mixture components or volatility path), its own \rho and \mu, and its own copy of prior. The series are coupled only through the multinomial likelihood, and each series is updated with the other categories held at their current values. A row with N_t = 0 carries no information about the shares and is handled like a missing observation. Zero inflation is not available for this family (zeros must be "none"). Running time grows linearly in K - 1.

Zero handling. Under zeros = "inflated" a latent gate decides, for each observed zero, whether it is structural (gate closed) or an ordinary sampling zero produced by the Poisson/binomial process (gate open); the gate-open probability is a single parameter that is constant over time. See structural_zero_prob(). Under zeros = "missing" the observed zeros are treated as missing values.

Forecasting. Forecasts are obtained by forward simulation from the posterior draws, so they can be computed after fitting for any horizon with forecast.dynamic_fit(). Setting horizon = H additionally simulates an H-step forecast right after sampling and stores it in the fit (under the same seed).

Reproducibility. With a seed, the sampler and the fit-time forecast run with that seed, and the previous state of the global random number generator is restored afterwards.

Value

An object of class "dynamic_fit": a list with three elements.

draws

Posterior draws only. For the Poisson and binomial families a list of matrices/vectors with one row (or element) per stored draw:

z

Latent states aligned with the observations (⁠draws x n⁠).

z0

The initial latent state, one period before the first observation, which carries the init_mean / init_var prior.

sig2

Variance of the increment leading into each z_t (⁠draws x n⁠); the first column is the increment from z0.

fitted, yrep

Unconditional fitted means and posterior predictive replicates; under zeros = "inflated" they include the zero-inflation gate.

fitted_open, yrep_open

Their conditional-on-gate-open counterparts: the latent-implied mean, and a replicate drawn straight from the observation model.

gate, pi_open

Zero-inflation gate indicators and gate-open probability; NULL unless zeros = "inflated".

rho, mu

AR(1) coefficient (1 under the random walk) and drift/intercept (0 unless include_mu = TRUE).

innov_var

A representative innovation variance whose definition depends on innovations: the estimated constant increment variance \sigma^2 for "gaussian"; the marginal increment variance for "t" (\sigma^2 \nu / (\nu - 2)) and "mixture" (\sigma^2 \sum_h \eta_h \sigma_h^2); and the average of the per-increment variances \exp(h_t) over the in-sample transitions for "sv". For "t", if a draw has \nu \le 2 (possible only when df_min <= 2), the undefined variance factor is replaced by the convention 3.

scale

The overall innovation scale \sigma^2 (NULL for "sv").

nu

Student-t degrees of freedom (NULL unless innovations = "t").

mix_weight, mix_var

Mixture weights and component variances (\sigma^2 \sigma_h^2), ⁠draws x mix_components⁠ (NULL unless innovations = "mixture"). The components are exchangeable and not identified individually. Within each draw they are stored in order of increasing variance, which is a labelling convention rather than an identification, so per-component summaries are meaningful only if the components are clearly separated. Label-invariant quantities, such as innov_var and the forecasts, are unaffected.

sv_mu, sv_phi, sv_sigma

Level, persistence and volatility of the AR(1) log-variance process (NULL unless innovations = "sv").

forecast_z, forecast_y

Latent and response forecasts when horizon >= 1 (NULL otherwise); forecast_y is unconditional, i.e. includes the gate.

Use yrep, not yrep_open, for posterior predictive checks; without zero inflation the conditional and unconditional pairs are identical. See predict.dynamic_fit().

For the multinomial family the same components carry an extra trailing category dimension, named by the category labels: latent quantities (z, sig2, forecast_z, mix_weight, mix_var) are arrays with one slice per non-baseline category, and z0, rho, mu, innov_var, scale, nu and the ⁠sv_*⁠ parameters are ⁠draws x (K - 1)⁠ matrices; response-scale quantities (fitted – the expected counts N_t p_{t,k} –, yrep, forecast_y, and the additional fitted_prob / forecast_prob holding the category shares) are ⁠draws x time x K⁠ arrays including the baseline, in the original column order. fitted_open, yrep_open, gate and pi_open are NULL.

data

The observed inputs: y, trials (the row totals for the multinomial family), the series length n, and the resolved per-observation offset. For the multinomial family also K, the categories (column labels), and the baseline index and baseline_name.

spec

The model and MCMC specification: family, innovations, latent_dynamics, include_mu, zeros, prior, nsave, nburn, thin, horizon, and the forecast-period inputs forecast_offset / forecast_trials used for the stored forecast (NULL when not applicable).

See Also

dynamic_prior(), forecast.dynamic_fit(), predict.dynamic_fit(), plot_fitted(), structural_zero_prob(), simulate_dynamic_multinomial()

Examples

sim <- simulate_dynamic_poisson(n = 60, sigma = 0.2, log_rate0 = 2, seed = 1)
fit <- fit_dynamic_model(sim$y, family = "poisson", nsave = 300, nburn = 200,
                         seed = 1)
summary(fit)
forecast(fit, horizon = 5)

# multinomial choice counts: three categories, baseline chosen automatically
simm <- simulate_dynamic_multinomial(n = 40, sigma = 0.15, trials = 100,
                                     alr0 = c(-1, -0.5), seed = 2)
fitm <- fit_dynamic_model(simm$y, family = "multinomial",
                          nsave = 200, nburn = 100, seed = 2)
fitm
summary(fitm)$params

DynCount documentation built on Sept. 28, 2026, 5:10 p.m.