The joint covariance of the fixed and random coefficients

The covariance of everything the fit estimates, beta, betad, theta and the random-effect coefficients b together, in one matrix. It is what a delta method over a fitted CURVE needs and what vcov() cannot return: a penalized smooth's wiggly part is a random-effect block even when the smooth is a population term, so a covariance that stops at the fixed effects covers none of it.

Usage

frm_joint_cov(object)

Arguments

object: A fitted frmtmb_fit.

Value

A list with

V: the p x p joint covariance.
names: length p; the PARAMETER COMPONENT each row belongs to ("beta", "betad", "b", "theta", ...), which is what frm_lp_basis()$coef_pos indexes.
labels: length p; one label per row, beta.<coefficient> for a fixed effect and b.<block>.<level> for a random one.
lost_pos: the rows of V that belong to a parameter whose standard error is lost (see summary()); integer(0) when none is.

On a fit whose fixed design is ill-conditioned, root and Vrest hold V as Vrest + root root'; see frm_lp_basis() for why a product must be formed from them. They are NULL otherwise. The other elements (null, units) are internal. A fit with no random effects has no joint precision, and V is then the fixed-effect covariance sdreport()$cov.fixed; names says so. A fit with a latent residual memoizes the sparse factor of its joint precision instead, and V is formed from it here: about 8 p^2 bytes, and about four times that while it is formed, refused above getOption("frmtmb.latres_dense_rows", 20000) rows.

A parameter whose standard error is lost has NaN in its row and column, as in vcov(). The rest of V propagates the uncertainty of the outer parameters along the directions the fit determines and none along the lost ones. When each lost parameter is a direction of its own, that is the covariance with the lost parameters held at their estimates. When the loss is one direction that mixes several parameters, a combination of them orthogonal to it keeps its variance: on a gr(g, by = f) fit that lost one direction over two sds, a slope effect's variance is 7.1e-4 here against 1.7e-6 with both sds held, and the prediction standard errors differ by at most 0.4 percent.

Details

vcov(full = TRUE) returns the OUTER parameter vector's covariance and its row names are documented to be exactly confint()'s, so b is not in it and cannot be added. frm_lp_basis() is the accessor for a design at newdata; this is the accessor for the covariance those coefficients have.

The result is memoized on the fit, so the joint-precision solve is paid once however many curves are drawn from it. It is also the ONLY route to the covariance of an AUTOSCALED fit (frmtmb_control(autoscale = TRUE)): a fresh RTMB::sdreport(getJointPrecision = TRUE) goes round the reparameterization and returns a covariance built on the unscaled Hessian, which is a different matrix and is not marked as one.

Examples

set.seed(1)
dd <- data.frame(x = rnorm(120), g = factor(rep(1:12, each = 10)))
dd$y <- rnorm(120, 1 + 2 * dd$x + rnorm(12, 0, 0.5)[dd$g], 0.4)
fit <- frm(bf(y ~ x + (1 | g)), data = dd)
jc <- frm_joint_cov(fit)
dim(jc$V)
#> [1] 16 16
table(jc$names)
#> 
#>     b  beta betad theta 
#>    12     2     1     1 
head(jc$labels)
#> [1] "beta.(Intercept)"        "beta.x"                 
#> [3] "betad.sigma_(Intercept)" "b.1 | g.1.(Intercept)"  
#> [5] "b.1 | g.2.(Intercept)"   "b.1 | g.3.(Intercept)"