Skip to contents

Extract a report-ready table of Sigma or correlation entries

Source: R/extract-sigma-table.R
extract_Sigma_table.Rd

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.

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". 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.

See also

extract_Sigma() for the underlying matrix extractor; extract_correlations() for pairwise correlation intervals; bootstrap_Sigma() for bootstrap uncertainty.

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`: the per-row diagonal at "site_species" 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 that
#>   diagonal (the folded Psi of `latent(..., unique = TRUE)`, or an explicit
#>   `indep()` / `unique()` at this grouping).
## 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.0747438    NA    NA            none            none latent     TRUE
#> 2   0.4725290    NA    NA            none            none latent    FALSE
#> 3   0.7320134    NA    NA            none            none latent     TRUE
#> 4  -0.4607182    NA    NA            none            none latent    FALSE
#> 5  -0.3126950    NA    NA            none            none latent    FALSE
#> 6   0.6305749    NA    NA            none            none latent     TRUE
#> 7   0.2279525    NA    NA            none            none latent    FALSE
#> 8   0.1547141    NA    NA            none            none latent    FALSE
#> 9  -0.1508470    NA    NA            none            none latent    FALSE
#> 10  0.5400930    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.5327414    NA    NA            none            none correlation    FALSE
#> 3   1.0000000    NA    NA            none            none correlation     TRUE
#> 4  -0.5596475    NA    NA            none            none correlation    FALSE
#> 5  -0.4602492    NA    NA            none            none correlation    FALSE
#> 6   1.0000000    NA    NA            none            none correlation     TRUE
#> 7   0.2991972    NA    NA            none            none correlation    FALSE
#> 8   0.2460572    NA    NA            none            none correlation    FALSE
#> 9  -0.2584843    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.69621429      0.3785295
#> 2 trait_2 0.32071108      0.4113023
#> 3 trait_3 0.30487924      0.3256956
#> 4 trait_4 0.07463553      0.4654574
## 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.4725290  0.14138034  0.75241558       bootstrap        provided latent
#> 2 -0.4607182 -0.70273318 -0.26710813       bootstrap        provided latent
#> 3 -0.3126950 -0.57235512 -0.08811883       bootstrap        provided latent
#> 4  0.2279525  0.01947117  0.48383929       bootstrap        provided latent
#> 5  0.1547141  0.01141523  0.29569109       bootstrap        provided latent
#> 6 -0.1508470 -0.32990832 -0.01355431       bootstrap        provided latent
#>   diagonal triangle
#> 1    FALSE    upper
#> 2    FALSE    upper
#> 3    FALSE    upper
#> 4    FALSE    upper
#> 5    FALSE    upper
#> 6    FALSE    upper
# }