Fit a varlasso model to multivariate time series data

fit is the primary function for fitting models with varlasso.

fit(
  data,
  shared_q = NULL,
  shared_r = NULL,
  fixed_r = NULL,
  est_trend = FALSE,
  varss = TRUE,
  b_diag = list(mu = 0.7, sd = 1),
  b_diag_min = 0,
  b_offdiag = list(mu = 0, sd = 1),
  off_diag_prior = c("normal", "student-t", "laplace", "hs", "normal_pp"),
  sigma_pro = list(mu = 0, sd = 1),
  sigma_obs = list(mu = 0, sd = 1),
  nu = NULL,
  sigma_scale = list(df = 3, sd = 2),
  hs = list(df = 1, df_global = 1, df_slab = 4, scale_global = 1, scale_slab = 2),
  unique_reg = FALSE,
  pars_to_monitor = NULL,
  iter = 1000,
  warmup = floor(iter/2),
  thin = 1,
  chains = 3,
  save_log_lik = FALSE,
  estimation = c("sampling", "optimizing", "vb"),
  est_hessian = FALSE,
  prior_only = FALSE,
  ...
)

Arguments

data

The data, as a long format dataframe

shared_q

Optional vector with integers indicating which process variance parameters are shared; defaults to constant process variances shared across species.

shared_r

Optional vector with integers indicating which observation variance parameters are shared; defaults to shared observation variances across species.

fixed_r

Optional scalar or vector of values to be passed in for a model where observation error variances are not estimated. Defaults to NULL, and them being estimated. Cannot be set to zero.

est_trend

Whether or not to estimate a species-specific trend, defaults to FALSE

varss

Boolean, whether to fit a VARSS model (defaults to true) or VAR model

b_diag

A 2 element list specifying the mean and sd of a normal prior on the diagonal elements. These elements are mu (default= 0.7) and sd (default = 1)

b_diag_min

The minimum value of the B diagonal, defaults to 0. Setting this lower, e.g. -1, allows for oscillations

b_offdiag

A 2 element list specifying the mean and sd of a normal prior on the off-diagonal elements. These elements are mu (default= 0) and sd (default = 1)

off_diag_prior

Character string denoting which prior to use for off-diagonal elements. Can be "normal" (independent priors estimated for each) parameter, "normal_pp" (normal prior with partial pooling, where model is hierarchical with an estimated standard deviation), "student-t" (hierarchical student-t model used with estimated scale and df), "laplace" (double-exponential model used with estimated scale) or "hs" (regularized horseshoe prior used, following parameterization in rstanarm)

sigma_pro

A 2 element list specifying the mean and sd of a half Student-t prior on the process variance standard deviation. Elements of the list are mu (default= 0) and sd (default = 1), and the prior df is fixed at 3

sigma_obs

A 2 element list specifying the mean and sd of a half Student-t prior on the process variance standard deviation. Elements of the list are mu (default= 0) and sd (default = 1), and the prior df is fixed at 3

nu

Optional known parameter (df) for Student-t priors

sigma_scale

A 2 element list specifying the df and sd parameters for a half Student-t prior on the scale parameter. This is only used for the Laplace and Student-t priors; default values are 3 for df and 2 for sd

hs

A list containing the hyperparameters of the horseshoe prior. Default values are df (df of the Student-t prior for local shrinkage, defaults to 1), df_global (Student-t df of the global shrinkage, defaults to 1), df_slab (df of the Student-t prior on the slab regularization, defaults to 4), scale_global (Scale for Student-t prior on global shrinkage, defaults to 1), scale_slab (Scale of the Student-t prior on shrikage parameter, defaults to 2).

unique_reg

Boolean, whether to estimate unique regularization parameters for the Student-t or Laplace models (defaults to false)

pars_to_monitor

Optional string vector of which parameters to monitor. By default this is NULL and only most important parameters are monitored (variances, B matrix, etc). This can be the word "all", in which case all parameters are saved (this may be good for loo::loo_moment_match() for example).

iter

Number of MCMC iterations, defaults to 1000

warmup

Warmup / burn in phase, defaults to 500

thin

MCMC thin, defaults to 1

chains

MCMC chains, defaults to 3

save_log_lik

Whether to save log_lik, defaults to FALSE for speed

estimation

Whether to do MCMC sampling ("sampling"; default), maximum posterior estimation ("optimizing") or Stan's variational algrotim ("vb")

est_hessian

Whether to estimate hessian matrix if doing MAP estimation, defaults to FALSE

prior_only

Whether to only generate samples from the prior distribution, defaults to FALSE

...

Extra arguments to pass to sampling

Value

an object of class 'stanfit'

Examples

# \donttest{

iter <- 50
warmup <- 20
thin <- 1
chains <- 1

n_ts <- 4
n_time <- 70
n_burn <- 50
xx <- matrix(0, n_time, n_ts)
set.seed(123)
xx[1, ] <- rnorm(n_ts)
B <- matrix(rnorm(n = n_ts * n_ts, mean = 0, sd = 0.1), n_ts, n_ts)
diag(B) <- runif(n_ts, min = 0.7, max = 0.9)

for (i in 2:n_time) {
  xx[i, ] <- B %*% xx[i - 1, ] + rnorm(n_ts, mean = 0, sd = 0.04)
}
xx <- xx[-c(1:n_burn), ]

yy <- xx + rnorm(n = nrow(xx) * ncol(xx), mean = 0, sd = 0.02)
df <- data.frame(
  "time" = rep(1:nrow(yy), ncol(yy)),
  "species" = sort(rep(1:ncol(yy), nrow(yy))),
  "y" = c(yy)
)

# fit basic model with normal independent priors
m <- fit(data = df, chains = chains, iter = iter, warmup = warmup, thin = thin)
#> 
#> SAMPLING FOR MODEL 'varlasso' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.0001 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: WARNING: There aren't enough warmup iterations to fit the
#> Chain 1:          three stages of adaptation as currently configured.
#> Chain 1:          Reducing each adaptation stage to 15%/75%/10% of
#> Chain 1:          the given number of warmup iterations:
#> Chain 1:            init_buffer = 3
#> Chain 1:            adapt_window = 15
#> Chain 1:            term_buffer = 2
#> Chain 1: 
#> Chain 1: Iteration:  1 / 50 [  2%]  (Warmup)
#> Chain 1: Iteration:  5 / 50 [ 10%]  (Warmup)
#> Chain 1: Iteration: 10 / 50 [ 20%]  (Warmup)
#> Chain 1: Iteration: 15 / 50 [ 30%]  (Warmup)
#> Chain 1: Iteration: 20 / 50 [ 40%]  (Warmup)
#> Chain 1: Iteration: 21 / 50 [ 42%]  (Sampling)
#> Chain 1: Iteration: 25 / 50 [ 50%]  (Sampling)
#> Chain 1: Iteration: 30 / 50 [ 60%]  (Sampling)
#> Chain 1: Iteration: 35 / 50 [ 70%]  (Sampling)
#> Chain 1: Iteration: 40 / 50 [ 80%]  (Sampling)
#> Chain 1: Iteration: 45 / 50 [ 90%]  (Sampling)
#> Chain 1: Iteration: 50 / 50 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.004481 seconds (Warm-up)
#> Chain 1:                0.612224 seconds (Sampling)
#> Chain 1:                0.616705 seconds (Total)
#> Chain 1: 
#> Warning: There were 16 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
#> https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
#> Warning: There were 1 chains where the estimated Bayesian Fraction of Missing Information was low. See
#> https://mc-stan.org/misc/warnings.html#bfmi-low
#> Warning: Examine the pairs() plot to diagnose sampling problems
#> Warning: The largest R-hat is 2.1, indicating chains have not mixed.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#r-hat
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#bulk-ess
#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#tail-ess
# fit same model with shrinkage using student-t priors
m_t <- fit(
  data = df, chains = chains, iter = iter, warmup = warmup, thin = thin,
  off_diag_prior = "student-t"
)
#> 
#> SAMPLING FOR MODEL 'varlasso' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 6.2e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.62 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: WARNING: There aren't enough warmup iterations to fit the
#> Chain 1:          three stages of adaptation as currently configured.
#> Chain 1:          Reducing each adaptation stage to 15%/75%/10% of
#> Chain 1:          the given number of warmup iterations:
#> Chain 1:            init_buffer = 3
#> Chain 1:            adapt_window = 15
#> Chain 1:            term_buffer = 2
#> Chain 1: 
#> Chain 1: Iteration:  1 / 50 [  2%]  (Warmup)
#> Chain 1: Iteration:  5 / 50 [ 10%]  (Warmup)
#> Chain 1: Iteration: 10 / 50 [ 20%]  (Warmup)
#> Chain 1: Iteration: 15 / 50 [ 30%]  (Warmup)
#> Chain 1: Iteration: 20 / 50 [ 40%]  (Warmup)
#> Chain 1: Iteration: 21 / 50 [ 42%]  (Sampling)
#> Chain 1: Iteration: 25 / 50 [ 50%]  (Sampling)
#> Chain 1: Iteration: 30 / 50 [ 60%]  (Sampling)
#> Chain 1: Iteration: 35 / 50 [ 70%]  (Sampling)
#> Chain 1: Iteration: 40 / 50 [ 80%]  (Sampling)
#> Chain 1: Iteration: 45 / 50 [ 90%]  (Sampling)
#> Chain 1: Iteration: 50 / 50 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.309595 seconds (Warm-up)
#> Chain 1:                0.597926 seconds (Sampling)
#> Chain 1:                0.907521 seconds (Total)
#> Chain 1: 
#> Warning: There were 6 divergent transitions after warmup. See
#> https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: There were 2 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
#> https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
#> Warning: Examine the pairs() plot to diagnose sampling problems
#> Warning: The largest R-hat is 1.53, indicating chains have not mixed.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#r-hat
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#bulk-ess
#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#tail-ess
# fit same model using student-t priors and known df
m_t_known <- fit(
  data = df, chains = chains, iter = iter, warmup = warmup, thin = thin,
  off_diag_prior = "student-t", nu = 4
)
#> 
#> SAMPLING FOR MODEL 'varlasso' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 6.2e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.62 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: WARNING: There aren't enough warmup iterations to fit the
#> Chain 1:          three stages of adaptation as currently configured.
#> Chain 1:          Reducing each adaptation stage to 15%/75%/10% of
#> Chain 1:          the given number of warmup iterations:
#> Chain 1:            init_buffer = 3
#> Chain 1:            adapt_window = 15
#> Chain 1:            term_buffer = 2
#> Chain 1: 
#> Chain 1: Iteration:  1 / 50 [  2%]  (Warmup)
#> Chain 1: Iteration:  5 / 50 [ 10%]  (Warmup)
#> Chain 1: Iteration: 10 / 50 [ 20%]  (Warmup)
#> Chain 1: Iteration: 15 / 50 [ 30%]  (Warmup)
#> Chain 1: Iteration: 20 / 50 [ 40%]  (Warmup)
#> Chain 1: Iteration: 21 / 50 [ 42%]  (Sampling)
#> Chain 1: Iteration: 25 / 50 [ 50%]  (Sampling)
#> Chain 1: Iteration: 30 / 50 [ 60%]  (Sampling)
#> Chain 1: Iteration: 35 / 50 [ 70%]  (Sampling)
#> Chain 1: Iteration: 40 / 50 [ 80%]  (Sampling)
#> Chain 1: Iteration: 45 / 50 [ 90%]  (Sampling)
#> Chain 1: Iteration: 50 / 50 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.00598 seconds (Warm-up)
#> Chain 1:                0.006516 seconds (Sampling)
#> Chain 1:                0.012496 seconds (Total)
#> Chain 1: 
#> Warning: There were 30 divergent transitions after warmup. See
#> https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: Examine the pairs() plot to diagnose sampling problems
# fit same model using laplace priors
m_lp <- fit(
  data = df, chains = chains, iter = iter, warmup = warmup, thin = thin,
  off_diag_prior = "laplace"
)
#> 
#> SAMPLING FOR MODEL 'varlasso' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 4.5e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.45 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: WARNING: There aren't enough warmup iterations to fit the
#> Chain 1:          three stages of adaptation as currently configured.
#> Chain 1:          Reducing each adaptation stage to 15%/75%/10% of
#> Chain 1:          the given number of warmup iterations:
#> Chain 1:            init_buffer = 3
#> Chain 1:            adapt_window = 15
#> Chain 1:            term_buffer = 2
#> Chain 1: 
#> Chain 1: Iteration:  1 / 50 [  2%]  (Warmup)
#> Chain 1: Iteration:  5 / 50 [ 10%]  (Warmup)
#> Chain 1: Iteration: 10 / 50 [ 20%]  (Warmup)
#> Chain 1: Iteration: 15 / 50 [ 30%]  (Warmup)
#> Chain 1: Iteration: 20 / 50 [ 40%]  (Warmup)
#> Chain 1: Iteration: 21 / 50 [ 42%]  (Sampling)
#> Chain 1: Iteration: 25 / 50 [ 50%]  (Sampling)
#> Chain 1: Iteration: 30 / 50 [ 60%]  (Sampling)
#> Chain 1: Iteration: 35 / 50 [ 70%]  (Sampling)
#> Chain 1: Iteration: 40 / 50 [ 80%]  (Sampling)
#> Chain 1: Iteration: 45 / 50 [ 90%]  (Sampling)
#> Chain 1: Iteration: 50 / 50 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.003542 seconds (Warm-up)
#> Chain 1:                0.00315 seconds (Sampling)
#> Chain 1:                0.006692 seconds (Total)
#> Chain 1: 
#> Warning: There were 30 divergent transitions after warmup. See
#> https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: Examine the pairs() plot to diagnose sampling problems
# fit same model using horseshoe priors
m_hs <- fit(
  data = df, chains = chains, iter = iter, warmup = warmup, thin = thin,
  off_diag_prior = "hs"
)
#> 
#> SAMPLING FOR MODEL 'varlasso' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 6e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.6 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: WARNING: There aren't enough warmup iterations to fit the
#> Chain 1:          three stages of adaptation as currently configured.
#> Chain 1:          Reducing each adaptation stage to 15%/75%/10% of
#> Chain 1:          the given number of warmup iterations:
#> Chain 1:            init_buffer = 3
#> Chain 1:            adapt_window = 15
#> Chain 1:            term_buffer = 2
#> Chain 1: 
#> Chain 1: Iteration:  1 / 50 [  2%]  (Warmup)
#> Chain 1: Iteration:  5 / 50 [ 10%]  (Warmup)
#> Chain 1: Iteration: 10 / 50 [ 20%]  (Warmup)
#> Chain 1: Iteration: 15 / 50 [ 30%]  (Warmup)
#> Chain 1: Iteration: 20 / 50 [ 40%]  (Warmup)
#> Chain 1: Iteration: 21 / 50 [ 42%]  (Sampling)
#> Chain 1: Iteration: 25 / 50 [ 50%]  (Sampling)
#> Chain 1: Iteration: 30 / 50 [ 60%]  (Sampling)
#> Chain 1: Iteration: 35 / 50 [ 70%]  (Sampling)
#> Chain 1: Iteration: 40 / 50 [ 80%]  (Sampling)
#> Chain 1: Iteration: 45 / 50 [ 90%]  (Sampling)
#> Chain 1: Iteration: 50 / 50 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.004482 seconds (Warm-up)
#> Chain 1:                0.720074 seconds (Sampling)
#> Chain 1:                0.724556 seconds (Total)
#> Chain 1: 
#> Warning: There were 17 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
#> https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
#> Warning: There were 1 chains where the estimated Bayesian Fraction of Missing Information was low. See
#> https://mc-stan.org/misc/warnings.html#bfmi-low
#> Warning: Examine the pairs() plot to diagnose sampling problems
#> Warning: The largest R-hat is 2.1, indicating chains have not mixed.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#r-hat
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#bulk-ess
#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#tail-ess

# include example for just sampling from priors
m <- fit(
  data = df, chains = chains, iter = iter,
  warmup = warmup, thin = thin, prior_only = TRUE
)
#> 
#> SAMPLING FOR MODEL 'varlasso' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 4.1e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.41 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: WARNING: There aren't enough warmup iterations to fit the
#> Chain 1:          three stages of adaptation as currently configured.
#> Chain 1:          Reducing each adaptation stage to 15%/75%/10% of
#> Chain 1:          the given number of warmup iterations:
#> Chain 1:            init_buffer = 3
#> Chain 1:            adapt_window = 15
#> Chain 1:            term_buffer = 2
#> Chain 1: 
#> Chain 1: Iteration:  1 / 50 [  2%]  (Warmup)
#> Chain 1: Iteration:  5 / 50 [ 10%]  (Warmup)
#> Chain 1: Iteration: 10 / 50 [ 20%]  (Warmup)
#> Chain 1: Iteration: 15 / 50 [ 30%]  (Warmup)
#> Chain 1: Iteration: 20 / 50 [ 40%]  (Warmup)
#> Chain 1: Iteration: 21 / 50 [ 42%]  (Sampling)
#> Chain 1: Iteration: 25 / 50 [ 50%]  (Sampling)
#> Chain 1: Iteration: 30 / 50 [ 60%]  (Sampling)
#> Chain 1: Iteration: 35 / 50 [ 70%]  (Sampling)
#> Chain 1: Iteration: 40 / 50 [ 80%]  (Sampling)
#> Chain 1: Iteration: 45 / 50 [ 90%]  (Sampling)
#> Chain 1: Iteration: 50 / 50 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.007243 seconds (Warm-up)
#> Chain 1:                0.008387 seconds (Sampling)
#> Chain 1:                0.01563 seconds (Total)
#> Chain 1: 
#> Warning: The largest R-hat is 1.31, indicating chains have not mixed.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#r-hat
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#bulk-ess
#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#tail-ess
# }