Extract experimental multivariate covariance standard errors — covariance_standard

Usage

covariance_standard_errors(object, ...)

Arguments

object: A fitted model object.
...: Reserved for future arguments.

Value

A named list of standard-error matrices genetic_covariance, residual_covariance, genetic_correlation, residual_correlation, and a per-trait heritability SE vector, for hsquared_fit objects that contain them.

Details

covariance_standard_errors() returns experimental large-sample (delta-method) standard errors for the multivariate genetic/residual covariance and correlation matrices and per-trait h^2, for an opt-in unstructured multivariate fit, when the engine returned them. The SEs come from the observed information: a central finite-difference Hessian of the multivariate REML log-likelihood with fd_step = 1e-4, not the univariate average-information matrix. See the information-matrix table on Can I fit and report this?.

Heavy caveats (engine multivariate REML validation row, partial): the strict per-seed recovery gate is still a non-pass (7/12 unstructured seeds in the updated study), but a 12-seed bias/MCSE study (twin HSquared.jl#78) shows no detectable bias – all six covariance parameters have |bias| <= 2*MCSE (largest 0.84*MCSE), a low-power non-rejection consistent with an unbiased estimator (not a proof), with EBV accuracy ~ 0.90 in both traits; a cold-start replication (#79) reaches the same optimum on all 12 seeds, so it is not a warm-start artifact. The per-seed gate failures reflect sampling variance of the estimated G at this design, not a detected bias. These SEs remain asymptotic, REML-only, unstructured-only (the engine refuses structured / factor-analytic fits, whose loadings are rotation-nonidentified), omitted at a flat/boundary optimum, not coverage-calibrated, with no external comparator, and not a validated capability.

The returned genetic_correlation SE is on the correlation (r) scale. R does not return a genetic-correlation interval. estimate +/- 1.96 * se is not an extractor and can leave (-1, 1) when |r_g| is large. Julia HSquared.genetic_correlation_interval() builds a Fisher-z interval (z = atanh(r), se_z = se_r / (1 - r^2)); that helper is not wired to R. When |r_g| is at +-1, the engine treats the optimum as flat or boundary and omits these SEs. print() can still say converged = TRUE.