Skip to contents

Multivariate split-R-hat for dynamic_beta coefficient paths

Source: R/rhat_dynamic.R
rhat_dynamic_beta.Rd

Given one or more multi-chain lame fits whose $BETA is [n_iter, p, T], computes the multivariate R-hat per coefficient \(k\) treating the length-\(T\) path as one multivariate observation per iteration.

Usage

rhat_dynamic_beta(fit_list, coefs = NULL)

Arguments

fit_list

A list of fitted lame objects from lame_parallel(..., chains = K), or any ame_chain_list produced by re-running lame() with different seeds.

coefs

Optional character vector of coefficient names to subset (matches dimnames(fit$BETA)[[2]]).

Value

Data frame with one row per coefficient: coef, rhat_mvt, rhat_max_univariate, n_chains, n_iter_per_chain, n_periods. rhat_mvt is the Brooks-Gelman multivariate statistic; rhat_max_univariate is the max over per-(k,t) split-R-hat values for comparison.

Details

For chain \(c\), let \(\beta^{(c)}_t \in \mathbb{R}^T\) be the path. With \(m\) chains and \(n\) iterations per chain, define $$W_k = \tfrac{1}{m}\sum_c S_c^{(k)}, \quad B_k = \tfrac{n}{m-1}\sum_c (\bar\beta_c^{(k)} - \bar\beta^{(k)})(\bar\beta_c^{(k)} - \bar\beta^{(k)})'$$ where \(S_c^{(k)}\) is the within-chain sample covariance of path \(k\) in chain \(c\). Then \(V_k = ((n-1)/n) W_k + ((m+1)/(mn)) B_k\) and \(\hat R_k^{mvt} = \sqrt{\lambda_{\max}(W_k^{-1} V_k)}\). For nearly degenerate covariances we add a tiny ridge to \(W_k\).

Examples

# \donttest{
data(YX_bin_list)
# note: pass seed = to lame() -- an external set.seed() does not vary
# the sampler, so it would produce identical chains
fit_list <- lapply(c(1L, 2L), function(s) {
  lame(YX_bin_list$Y, YX_bin_list$X, family = "binary", R = 0,
       dynamic_beta = "dyad", seed = s,
       nscan = 60, burn = 15, odens = 5, verbose = FALSE)
})
#> Warning: `family` = "binary" but `Y` contains values other than 0/1.
#> ℹ `Y` will be thresholded to `1 * (Y > 0)`; if you meant counts, use "poisson",
#>   or "ordinal"/"normal" as appropriate.
#> Warning: `family` = "binary" but `Y` contains values other than 0/1.
#> ℹ `Y` will be thresholded to `1 * (Y > 0)`; if you meant counts, use "poisson",
#>   or "ordinal"/"normal" as appropriate.
rhat_dynamic_beta(fit_list)
#>        coef  rhat_mvt rhat_max_univariate n_chains n_iter_per_chain n_periods
#> 1 intercept 0.9701234            1.197511        2               12         4
#> 2   X1_dyad 1.3518445            1.848351        2               12         4
#> 3   X2_dyad 2.9315631            1.991456        2               12         4
#> 4   X3_dyad 2.4477649            2.033285        2               12         4
# }