Skip to contents

Heritability, ICC, and repeatability accessors

Source: R/heritability.R
heritability.Rd

heritability(), icc(), and repeatability() are derived-quantity accessors for structured-Gaussian drmTMB fits. All three report the share of total variance carried by one structured random-effect component relative to a denominator on the working (log-SD) scale:

Usage

heritability(object, ...)

# S3 method for class 'drmTMB'
heritability(
  object,
  component = NULL,
  level = 0.95,
  method = c("delta", "profile"),
  ...
)

icc(object, ...)

# S3 method for class 'drmTMB'
icc(
  object,
  component = NULL,
  level = 0.95,
  method = c("delta", "profile"),
  ...
)

repeatability(object, ...)

# S3 method for class 'drmTMB'
repeatability(
  object,
  component = NULL,
  level = 0.95,
  method = c("delta", "profile"),
  ...
)

# S3 method for class 'drmTMB_julia'
heritability(object, ...)

# S3 method for class 'drmTMB_julia'
icc(object, ...)

# S3 method for class 'drmTMB_julia'
repeatability(object, ...)

Arguments

object

A drmTMB fit.

...

Reserved for future extractor options.

component

Optional character string naming the structured mean component (matching a name in object$sdpars$mu, e.g. "(1 | id)" or "phylo(1 | species)"; structured markers drop their non-grouping arguments such as tree =/Ainv =/K = from the label – see the @examples phylo case below). Required when the fit has more than one structured component; ignored (and unnecessary) when the fit has exactly one.

level

Confidence level for the Wald interval. Default 0.95.

method

Either "delta" (default) or "profile". "profile" is not implemented and always aborts.

Value

A one-row data frame (class drm_heritability) with columns quantity, component, estimate, se, lower, upper, level, and method.

Details

  • heritability(): sigma^2_focal / (sum of ALL structured mean-component variances + residual variance) – the comparative-biology "phylogenetic signal" definition.

  • icc() / repeatability(): sigma^2_focal / (sigma^2_focal + residual variance) – the classic focal-vs-residual intraclass correlation. repeatability() is an alias for icc().

heritability() and icc() return the same value only when the fit has a single structured mean component (no other components to net out of the denominator).

summary()'s derived rows, named total_variance_share and phylo_total_variance_share, use the total-variance (heritability-style) denominator, so they can differ from icc()/repeatability() here whenever a fit has two or more structured components. Their names make the denominator explicit so that one label does not refer to two different quantities.

Fits must be Gaussian, have a constant residual scale (sigma ~ 1) or a closed-form marginal residual variance, and have at least one structured mean random-effect component ((1 | group), phylo(...), animal(...), relmat(...), or spatial(...)); a component modelled with sd(group) ~ ... (location-scale-scale) does not define a single variance component and is refused. When a fit has more than one structured component, component must name one of them (see object$sdpars$mu for the available labels). When sigma additionally carries an ordinary random intercept, or a phylogenetic random intercept with a unit-diagonal correlation (the default phylo(...) route), the residual entry of the denominator is the marginal residual variance exp(2*b0 + 2*sum_k(omega_k^2)) – E[sigma^2], not the squared median exp(2*b0) – and its delta-method gradient is extended over those omega_k positions too. A random slope on sigma, a structured sigma effect without a verified unit-diagonal correlation (measured on the rows the design uses), or a fit that carries a sigma random effect while the log(sigma) soft clamp bent the assembled predictor for at least one observation (clamp_limited), has no closed-form marginal residual variance and is refused with a message naming why.

Standard errors and confidence intervals use a delta method on the working (log-SD) scale: a numeric gradient of the ratio in the log-SD parameters, combined with the fixed-parameter covariance (object$sdr$cov.fixed) via the quadratic form, then a Wald interval at level, clamped to [0, 1]. This is a delta-method approximation; not a coverage claim.

method = "profile" is accepted but is not implemented; it aborts naming method = "delta" as the supported option.

Examples

set.seed(20260525)
n_groups <- 20
n_per <- 8
grp <- factor(rep(seq_len(n_groups), each = n_per))
sd_g <- 1
sd_e <- 0.6
b_g <- rnorm(n_groups, sd = sd_g)
dat <- data.frame(y = 2 + b_g[grp] + rnorm(length(grp), sd = sd_e), grp = grp)
fit <- drmTMB(bf(y ~ 1 + (1 | grp), sigma ~ 1), data = dat)
icc(fit)
#> <icc> component = "(1 | grp)"
#> estimate: 0.7417
#> se: 0.06758 95% CI: [0.6093, 0.8742]
repeatability(fit)
#> <repeatability> component = "(1 | grp)"
#> estimate: 0.7417
#> se: 0.06758 95% CI: [0.6093, 0.8742]
heritability(fit)
#> <heritability> component = "(1 | grp)"
#> estimate: 0.7417
#> se: 0.06758 95% CI: [0.6093, 0.8742]

# \donttest{
if (requireNamespace("ape", quietly = TRUE)) {
  # A structured (phylo) component: its sdpars$mu / component label drops
  # the tree argument entirely and is just "phylo(1 | species)".
  set.seed(20260601)
  n_tip <- 20
  phy <- ape::rcoal(n_tip)
  phy$tip.label <- paste0("sp_", seq_len(n_tip))
  A <- ape::vcv(phy, corr = TRUE)
  u <- as.vector(t(chol(A)) %*% rnorm(n_tip)) * 0.9
  species <- factor(rep(phy$tip.label, each = 4), levels = phy$tip.label)
  phylo_dat <- data.frame(
    y = 2 + u[rep(seq_len(n_tip), each = 4)] +
      rnorm(length(species), sd = 0.6),
    species = species
  )
  phylo_fit <- drmTMB(
    bf(y ~ 1 + phylo(1 | species, tree = phy), sigma ~ 1),
    data = phylo_dat
  )
  icc(phylo_fit, component = "phylo(1 | species)")
}
#> <icc> component = "phylo(1 | species)"
#> estimate: 0.4374
#> se: 0.2074 95% CI: [0.03085, 0.8439]
# }