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, ported term-for-term from DRM.jl's src/heritability.jl (design doc docs/design/259-heritability-icc-repeatability.md). 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 accepted for signature parity only and always aborts in this slice.

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. Those rows were called repeatability and phylogenetic_signal until 2026-09-03; they were renamed so that one word no longer names two different quantities (D-213).

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; see docs/design/275-repeatability-scale-and-residual-variance.md.

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. See the design doc for the sanity-check evidence this slice carries.

method = "profile" is accepted for call-site parity with DRM.jl but is not implemented in this slice; 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]
# }