Extract a report-ready table of Sigma or correlation entries

extract_Sigma() returns matrices because that is the most convenient interactive representation. extract_Sigma_table() returns the same covariance target as one row per matrix entry so articles, tables, and plot helpers can work without hand-indexing matrices. It can also turn a bootstrap_Sigma() result into the same row schema with bootstrap interval columns filled in.

Usage

extract_Sigma_table(
  fit,
  level = "unit",
  part = c("total", "shared", "unique"),
  measure = c("covariance", "correlation"),
  entries = c("unique", "all", "upper", "lower", "offdiag", "diag"),
  link_residual = c("auto", "none")
)

Arguments

fit: A fit returned by gllvmTMB(), an admitted engine = "julia" bridge fit, or a bootstrap_Sigma() result.
level: Character vector of covariance levels, or "all" for every level present in the fit. Canonical levels are "unit", "unit_obs", "cluster", "cluster2", "phy", and "spatial"; named kernel_*() tiers are accepted by their fitted name; legacy aliases "B", "W", and "spde" are accepted. For gllvmTMB_julia objects, only "unit" is currently routed; "all" maps to that tier.
part: One of "total" (default), "shared", or "unique", passed to extract_Sigma().
measure: One of "covariance" (default) or "correlation". Correlation tables are available only for part = "total" because the report-ready correlation target is based on the full Lambda Lambda^T + Psi covariance.
entries: Which symmetric matrix entries to return. "unique" (default) returns the diagonal plus upper triangle, one row per unique estimand. "all" returns every cell, useful for heatmaps. "upper", "lower", "offdiag", and "diag" return the corresponding subsets. Ignored when part = "unique": that part always returns the T diagonal Psi rows, whatever entries is.
link_residual: Passed to extract_Sigma(). "auto" (default) adds family/link implicit residual variances to non-Gaussian trait diagonals; "none" returns only the fitted model-implied covariance.

Value

A data frame with one row per requested entry and stable columns: estimand, trait_i, trait_j, integer indices i and j, level, component, matrix, estimate, lower, upper, interval_method, interval_status, scale, diagonal, triangle, and validation_row (a short label for the validation scope of the extractor behind the table). For a fit, the interval columns lower and upper are NA with interval_method = "none" because this helper is point-estimate only. For a bootstrap_Sigma() input they are filled, with interval_method = "bootstrap" and interval_status = "provided".

Details

The helper reshapes the same covariance and correlation matrices returned by extract_Sigma(). A bootstrap_Sigma() object can supply interval columns. Julia bridge fits currently expose only point estimates for the ordinary unit tier. The helper does not create profile or Wald intervals.

The table is a point-estimate view over extract_Sigma(). It does not compute confidence intervals from a fitted model directly. Use extract_correlations() when you need pairwise correlation intervals, or bootstrap_Sigma() followed by extract_Sigma_table() when you need bootstrap intervals for Sigma or correlation matrix entries.

Examples

## Simulate 4 traits recorded on a few species at each of 50 sites, with
## correlated trait variation at two tiers: between sites (the "unit"
## tier) and between species records within a site (the "unit_obs"
## tier). A zero beta matrix turns off the simulator's environmental
## effects, which this model does not include.
set.seed(3)
sim <- simulate_site_trait(
  n_sites = 50, n_species = 10, n_traits = 4,
  mean_species_per_site = 3,
  beta = matrix(0, 4, 2),
  Lambda_B = matrix(c(0.8, 0.6, -0.5, 0.4), 4, 1),
  psi_B = rep(0.3, 4),
  Lambda_W = matrix(c(0.5, -0.4, 0.4, 0.3), 4, 1),
  sigma2_eps = 0.3
)
## One latent factor plus trait-specific variances at each tier. Each
## species record has one value per trait, so the within-site
## trait-specific variances also absorb the measurement noise
## (gllvmTMB prints a note saying so).
fit <- gllvmTMB(
  value ~ 0 + trait +
    latent(0 + trait | site, d = 1, unique = TRUE) +
    latent(0 + trait | site_species, d = 1, unique = TRUE),
  data     = sim$data,
  trait    = "trait",
  unit     = "site",
  unit_obs = "site_species"
)
#> ℹ Auto-suppressing `sigma_eps`: `indep(0 + trait | site_species)` is at the
#>   per-row level, so it already absorbs the observation residual.
#> • Fixed at 0.00122 (~1/1000 of sd(y)) to keep the Gaussian density
#>   well-defined; the row-level residual variance is fully captured by the
#>   per-row diagonal term.
## Between-site covariance: one row per unique entry (diagonal plus
## upper triangle), point estimates only
extract_Sigma_table(fit, level = "unit")
#>                       estimand trait_i trait_j i j level component matrix
#> 1  Sigma_unit[trait_1,trait_1] trait_1 trait_1 1 1  unit     total  Sigma
#> 2  Sigma_unit[trait_1,trait_2] trait_1 trait_2 1 2  unit     total  Sigma
#> 3  Sigma_unit[trait_2,trait_2] trait_2 trait_2 2 2  unit     total  Sigma
#> 4  Sigma_unit[trait_1,trait_3] trait_1 trait_3 1 3  unit     total  Sigma
#> 5  Sigma_unit[trait_2,trait_3] trait_2 trait_3 2 3  unit     total  Sigma
#> 6  Sigma_unit[trait_3,trait_3] trait_3 trait_3 3 3  unit     total  Sigma
#> 7  Sigma_unit[trait_1,trait_4] trait_1 trait_4 1 4  unit     total  Sigma
#> 8  Sigma_unit[trait_2,trait_4] trait_2 trait_4 2 4  unit     total  Sigma
#> 9  Sigma_unit[trait_3,trait_4] trait_3 trait_4 3 4  unit     total  Sigma
#> 10 Sigma_unit[trait_4,trait_4] trait_4 trait_4 4 4  unit     total  Sigma
#>      estimate lower upper interval_method interval_status  scale diagonal
#> 1   1.0747672    NA    NA            none            none latent     TRUE
#> 2   0.4725395    NA    NA            none            none latent    FALSE
#> 3   0.7320213    NA    NA            none            none latent     TRUE
#> 4  -0.4607257    NA    NA            none            none latent    FALSE
#> 5  -0.3127059    NA    NA            none            none latent    FALSE
#> 6   0.6305812    NA    NA            none            none latent     TRUE
#> 7   0.2279548    NA    NA            none            none latent    FALSE
#> 8   0.1547185    NA    NA            none            none latent    FALSE
#> 9  -0.1508505    NA    NA            none            none latent    FALSE
#> 10  0.5400961    NA    NA            none            none latent     TRUE
#>    triangle
#> 1  diagonal
#> 2     upper
#> 3  diagonal
#> 4     upper
#> 5     upper
#> 6  diagonal
#> 7     upper
#> 8     upper
#> 9     upper
#> 10 diagonal
## The same entries on the correlation scale
extract_Sigma_table(fit, level = "unit", measure = "correlation")
#>                   estimand trait_i trait_j i j level component matrix
#> 1  R_unit[trait_1,trait_1] trait_1 trait_1 1 1  unit     total      R
#> 2  R_unit[trait_1,trait_2] trait_1 trait_2 1 2  unit     total      R
#> 3  R_unit[trait_2,trait_2] trait_2 trait_2 2 2  unit     total      R
#> 4  R_unit[trait_1,trait_3] trait_1 trait_3 1 3  unit     total      R
#> 5  R_unit[trait_2,trait_3] trait_2 trait_3 2 3  unit     total      R
#> 6  R_unit[trait_3,trait_3] trait_3 trait_3 3 3  unit     total      R
#> 7  R_unit[trait_1,trait_4] trait_1 trait_4 1 4  unit     total      R
#> 8  R_unit[trait_2,trait_4] trait_2 trait_4 2 4  unit     total      R
#> 9  R_unit[trait_3,trait_4] trait_3 trait_4 3 4  unit     total      R
#> 10 R_unit[trait_4,trait_4] trait_4 trait_4 4 4  unit     total      R
#>      estimate lower upper interval_method interval_status       scale diagonal
#> 1   1.0000000    NA    NA            none            none correlation     TRUE
#> 2   0.5327446    NA    NA            none            none correlation    FALSE
#> 3   1.0000000    NA    NA            none            none correlation     TRUE
#> 4  -0.5596477    NA    NA            none            none correlation    FALSE
#> 5  -0.4602605    NA    NA            none            none correlation    FALSE
#> 6   1.0000000    NA    NA            none            none correlation     TRUE
#> 7   0.2991962    NA    NA            none            none correlation    FALSE
#> 8   0.2460622    NA    NA            none            none correlation    FALSE
#> 9  -0.2584881    NA    NA            none            none correlation    FALSE
#> 10  1.0000000    NA    NA            none            none correlation     TRUE
#>    triangle
#> 1  diagonal
#> 2     upper
#> 3  diagonal
#> 4     upper
#> 5     upper
#> 6  diagonal
#> 7     upper
#> 8     upper
#> 9     upper
#> 10 diagonal
## Split each trait's between-site variance into the part shared
## through the latent factor and the trait-specific part; the two add
## up to the diagonal of the first table
shared <- extract_Sigma_table(fit, level = "unit", part = "shared",
                              entries = "diag")
own <- extract_Sigma_table(fit, level = "unit", part = "unique",
                           entries = "diag")
data.frame(trait = shared$trait_i, shared = shared$estimate,
           trait_specific = own$estimate)
#>     trait     shared trait_specific
#> 1 trait_1 0.69621687      0.3785504
#> 2 trait_2 0.32072421      0.4112971
#> 3 trait_3 0.30488801      0.3256932
#> 4 trait_4 0.07463678      0.4654593
## Every cell of both tiers in one long table, handy for a heatmap:
## 2 tiers x 4 x 4 traits = 32 rows
all_cells <- extract_Sigma_table(fit, level = c("unit", "unit_obs"),
                                 entries = "all")
table(all_cells$level)
#> 
#>     unit unit_obs 
#>       16       16 
# \donttest{
## Bootstrap intervals: bootstrap_Sigma() refits the model to 50 data
## sets simulated from the fit, so this part takes tens of seconds.
## A small n_boot keeps the run time down and makes bootstrap_Sigma()
## warn that so few resamples give noisy bounds; real analyses need
## several hundred (the default is 999). These percentile intervals are
## not guaranteed to reach their nominal coverage.
boot <- bootstrap_Sigma(fit, n_boot = 50, level = "unit",
                        what = "Sigma", progress = FALSE)
#> Warning: `n_boot` = 50 gives noisy percentile bounds.
#> ℹ Percentile CIs are first-order accurate, so more draws buy precision, not
#>   correctness -- but below ~999 the Monte Carlo error in the endpoints is
#>   itself material.
#> ℹ `n_cores` >= 2 parallelises the refits; the cost here is wall-clock, not
#>   accuracy.
#> ℹ For Sigma diagonals, `profile_ci_total_variance()` needs no resampling at
#>   all.
extract_Sigma_table(boot, level = "unit", entries = "upper")
#>                      estimand trait_i trait_j i j level component matrix
#> 1 Sigma_unit[trait_1,trait_2] trait_1 trait_2 1 2  unit     total  Sigma
#> 2 Sigma_unit[trait_1,trait_3] trait_1 trait_3 1 3  unit     total  Sigma
#> 3 Sigma_unit[trait_2,trait_3] trait_2 trait_3 2 3  unit     total  Sigma
#> 4 Sigma_unit[trait_1,trait_4] trait_1 trait_4 1 4  unit     total  Sigma
#> 5 Sigma_unit[trait_2,trait_4] trait_2 trait_4 2 4  unit     total  Sigma
#> 6 Sigma_unit[trait_3,trait_4] trait_3 trait_4 3 4  unit     total  Sigma
#>     estimate       lower       upper interval_method interval_status  scale
#> 1  0.4725395  0.14138267  0.75242868       bootstrap        provided latent
#> 2 -0.4607257 -0.70274744 -0.26711826       bootstrap        provided latent
#> 3 -0.3127059 -0.57236175 -0.08812277       bootstrap        provided latent
#> 4  0.2279548  0.01947245  0.48384589       bootstrap        provided latent
#> 5  0.1547185  0.01141635  0.29570265       bootstrap        provided latent
#> 6 -0.1508505 -0.32991451 -0.01355529       bootstrap        provided latent
#>   diagonal triangle
#> 1    FALSE    upper
#> 2    FALSE    upper
#> 3    FALSE    upper
#> 4    FALSE    upper
#> 5    FALSE    upper
#> 6    FALSE    upper
# }