gllvmTMB() computes a TMB TMB::sdreport() at fitting time unless it is
told not to (control = gllvmTMBcontrol(se = FALSE)). Skipping it is
substantially faster, but it used to be a one-way door: the only way to get
standard errors afterwards was to fit the model again. standard_errors()
removes that door by computing the same sdreport() on demand from the
fitted object.
Arguments
- fit
A fitted multivariate model returned by
gllvmTMB().
Value
The same fit, with its sd_report field populated (and
sdreport_error cleared). Standard-error consumers such as
summary.gllvmTMB_multi(), getLV(), getREsd(), and
confint(method = "wald") work on the returned object.
R semantics matter here: the fit is returned, not modified in place.
Assign the result – fit <- standard_errors(fit) – or the standard
errors are discarded.
If the fit already has an sd_report, it is returned unchanged; the
calculation is not repeated.
Details
Fits made with non-unit likelihood weights use a weighted objective rather
than an ordinary likelihood. standard_errors() refuses those fits because
the package does not yet provide sandwich-calibrated uncertainty for them;
point estimates remain available from the fitted model.
When this helps
The saving is real but bounded: it is worth using when a workflow fits many
models and needs standard errors for only some of them (model selection,
simulation loops, a grid of candidate structures). For a single fit whose
summary you intend to read, the default se = TRUE is simpler and costs
nothing extra.
Same-session only
A TMB ADFun object holds external pointers into compiled memory. Those
pointers do not survive saveRDS() / readRDS() or a new R session, so
standard_errors() can only be called on a fit made in the current
session. This is a pre-existing property of every function in this package
that reuses fit$tmb_obj (profiling, bootstrapping, getLV()), not a new
restriction – but this function names it with a clear error rather than
failing obscurely. If you need standard errors from a saved fit, refit it
with se = TRUE.
What these standard errors are
Exactly the ones the fit would have produced at fitting time – the same
TMB sdreport() call, on the same converged parameter vector, to the bit.
Deferring the calculation does not change it, and does not change any
honesty caveat that already applies to Wald standard errors from this
package.
With finite loading ridge, the calculation includes the Gaussian loading prior precision at the fitted MAP point. These are local approximate posterior standard deviations conditional on the chosen ridge scale, model, and latent-axis orientation, not calibrated sampling standard errors.
See also
gllvmTMBcontrol() for the se argument that defers the
calculation; getREsd() and getLV() for accessors that read the
resulting sd_report.
Examples
## Simulate 40 sites, each visited twice. Four traits share one latent
## site axis, and each trait also has its own site-level variation.
## The repeat visits let the model separate that trait-specific site
## variation from visit-to-visit noise.
set.seed(2)
n_site <- 40
dat <- data.frame(site = factor(rep(seq_len(n_site), each = 2)))
z <- rnorm(n_site) # latent site scores
lambda <- c(0.9, 0.7, 0.6, -0.5) # loadings on that axis
for (j in 1:4) {
site_eff <- lambda[j] * z + rnorm(n_site, sd = 0.6)
dat[[paste0("y", j)]] <- site_eff[dat$site] +
rnorm(nrow(dat), sd = 0.5)
}
## Fit fast, decide later whether standard errors are needed.
fit <- gllvmTMB(
traits(y1, y2, y3, y4) ~ 1 + latent(1 | site, d = 1, unique = TRUE),
data = dat, family = gaussian(), unit = "site",
control = gllvmTMBcontrol(se = FALSE)
)
## Assign the result: standard_errors() returns an updated copy of the
## fit rather than changing `fit` in place.
fit <- standard_errors(fit)
## The Std.Err column of the fixed effects is now filled in (it was NA
## before standard_errors() ran).
summary(fit)
#> Stacked-trait gllvmTMB summary
#> Traits = 4, site = 40
#> Covstructs: latent_unit, indep_unit
#> ML log L = -380.093 convergence = 0 engine = Laplace
#>
#> Fixed effects:
#> Estimate Std.Err
#> traity1 0.104 0.195
#> traity2 -0.008 0.147
#> traity3 0.339 0.136
#> traity4 -0.161 0.133
#>
#> Between-unit trait correlation (R_B):
#> y1 y2 y3 y4
#> y1 1.000 0.455 0.527 -0.464
#> y2 0.455 1.000 0.586 -0.516
#> y3 0.527 0.586 1.000 -0.599
#> y4 -0.464 -0.516 -0.599 1.000
#>
#> Per-trait variance summaries:
#> comm_B
#> y1 0.409
#> y2 0.505
#> y3 0.680
#> y4 0.527
#>
#> For more, see: extract_Sigma(), extract_communality(),
#> extract_phylo_signal(), extract_proportions(), getLoadings(),
#> bootstrap_Sigma(), gllvmTMB_diagnose(), or plot(fit, type = ...).