Heritability tables with R-itable and clerkR

Note: This vignette requires the Ritable package (remotes::install_github("circadia-bio/R-itable")). Code is shown below but not executed in this build.

Overview

clerkR::tbl_heritability() is designed to accept the output of R-itable::herit_batch() directly. Column-name defaults match herit_batch() output exactly:

`tbl_heritability()` argument	Default	`herit_batch()` column
`metric`	`"trait"`	trait name
`h2`	`"h2"`	heritability estimate
`ci_low`	`"ci_lo"`	lower CI bound
`ci_high`	`"ci_hi"`	upper CI bound
`p`	`"pval"`	one-sided LRT p-value

The covariate model column is "covariates" in herit_batch() output — pass model = "covariates" to include it.

Setup

Install both packages from GitHub if you haven’t already:

# install.packages("remotes")
remotes::install_github("circadia-bio/R-itable")
remotes::install_github("circadia-bio/clerkR")

library(Ritable)
library(clerkR)

1. Simulate a family cohort

We simulate 80 nuclear families (2 parents + 3 offspring each) with three quantitative traits at different true h² values — BMI (~0.60), systolic blood pressure (~0.35), and HDL cholesterol (~0.50).

set.seed(2026)
n_fam <- 80

ped <- do.call(rbind, lapply(seq_len(n_fam), function(f) {
  base <- (f - 1L) * 5L
  data.frame(
    id  = base + 1:5,
    pat = c(0L, 0L, base + 1L, base + 1L, base + 1L),
    mom = c(0L, 0L, base + 2L, base + 2L, base + 2L),
    sex = c(1L, 2L, sample(1:2, 3, replace = TRUE))
  )
}))

simulate_trait <- function(ped, n_fam, h2, seed) {
  set.seed(seed)
  genetic  <- rep(rnorm(n_fam, 0, sqrt(h2)),      each = 5)
  residual <- rnorm(nrow(ped),         0, sqrt(1 - h2))
  genetic + residual
}

pheno_all <- data.frame(
  bmi = simulate_trait(ped, n_fam, h2 = 0.60, seed = 1),
  sbp = simulate_trait(ped, n_fam, h2 = 0.35, seed = 2),
  hdl = simulate_trait(ped, n_fam, h2 = 0.50, seed = 3)
)

off_rows <- ped$pat != 0L
dat <- data.frame(
  IID     = ped$id[off_rows],
  age     = round(runif(sum(off_rows), 20, 75)),
  sex_num = ped$sex[off_rows],
  bmi     = pheno_all$bmi[off_rows],
  sbp     = pheno_all$sbp[off_rows],
  hdl     = pheno_all$hdl[off_rows]
)
dat$age2 <- dat$age^2

2. Build GRM and run batch estimation

A <- build_grm(ped, study_ids = dat$IID)

res <- herit_batch(
  traits    = c("bmi", "sbp", "hdl"),
  grm       = A,
  data      = dat,
  covs_list = list(
    unadj = NULL,
    cov1  = c("age", "sex_num"),
    cov2  = c("age", "sex_num", "age2")
  )
)

# Inspect available covariate model labels
unique(res$covariates)
res

3. Render the full table

Pass column names as strings — defaults already match herit_batch() output, so only model, sigma2_a, and sigma2_e need to be specified:

res |>
  tbl_heritability(
    model    = "covariates",
    sigma2_a = "sigma2_a",
    sigma2_e = "sigma2_e",
    domains  = list("Cardiometabolic" = c("bmi", "sbp", "hdl")),
    fdr      = TRUE,
    output   = "gt"
  ) |>
  clerk_render(
    title    = "Heritability estimates",
    subtitle = "Narrow-sense h\u00b2 \u00b1 95% profile-likelihood CI",
    footnote = "FDR correction applied within each covariate model (BH). p-values are one-sided LRT with chi-squared(1) boundary correction."
  )

Filtering to a single model

adj_model <- unique(res$covariates)[length(unique(res$covariates))]

res[res$covariates == adj_model, ] |>
  tbl_heritability(
    sigma2_a = "sigma2_a",
    sigma2_e = "sigma2_e",
    fdr      = TRUE,
    output   = "gt"
  ) |>
  clerk_render(
    title    = paste0("Heritability estimates (", adj_model, ")"),
    footnote = "FDR correction applied across all traits (BH)."
  )

HTML output

res[res$covariates == adj_model, ] |>
  tbl_heritability(
    sigma2_a = "sigma2_a",
    sigma2_e = "sigma2_e",
    output   = "html"
  ) |>
  clerk_render()