Skip to contents

A vs C identifiability check (ACE model comparison)

Source: R/herit_ace.R
herit_ace.Rd

Fits four nested variance-component models for a single quantitative trait – sporadic (no familial effect), polygenic/AE (additive genetic only), household/CE (common environment only), and full ACE (A and C estimated jointly) – and returns their log-likelihoods together with the two likelihood-ratio tests needed to judge whether A and C are separately identifiable in the sample: CE vs sporadic (is there a familial signal at all?) and ACE vs AE (does C add anything once A is in the model?). This mirrors SOLAR Eclipse's house / polygenic -screen model-comparison workflow (see solar_full_analysis.tcl, Part 3).

Usage

herit_ace(
  trait,
  grm,
  household,
  data,
  covs = NULL,
  screen = TRUE,
  prob_level = 0.1,
  id_col = "IID",
  transform = TRUE,
  min_n = 80L,
  verbose = TRUE
)

Arguments

trait

Character string: name of the trait column in data.

grm

Numeric matrix: additive genetic relationship matrix, as returned by build_grm().

household

Numeric matrix: household indicator matrix, as returned by build_household(). Must have the same row/column names as grm.

data

Data frame containing id_col, trait, and any covs.

covs

Optional character vector of covariate column names in data (e.g. c("age", "sex", "age_sex", "age2", "age2_sex") for SOLAR's standard age/sex/agesex/age^2/age^2sex set – construct interaction and polynomial terms as ordinary columns before calling, same convention as herit_vc()). Covariates are standardised internally. Default NULL (no covariates, matching prior behaviour exactly).

screen

Logical. If covs is given, apply SOLAR-style automatic covariate screening (see Details) rather than including all of covs unconditionally. Default TRUE. Ignored if covs is NULL.

prob_level

Screening significance threshold. Default 0.10, matching SOLAR Eclipse's own default (confirmed from source, lib/solar.tcl proc polygenic: set probability_level 0.1).

id_col

Name of the individual ID column in data. Default "IID".

transform

Logical. Apply int_transform() to trait before fitting. Default TRUE, matching SOLAR Eclipse convention for skewed phenotypes.

min_n

Minimum number of complete observations required. Default 80L.

verbose

Logical. Print progress via cli. Default TRUE.

Value

A named list with elements:

trait

Trait name.

n

Sample size after dropping missing values.

loglik_sporadic, loglik_ae, loglik_ce, loglik_ace

Log-likelihoods of the four nested models (all fit with the same final covariate set – see Details).

h2_ae

Heritability under the AE model.

c2_ce

Common-environment proportion under the CE model.

h2_ace, c2_ace

Joint estimates under the full ACE model.

chisq_c_vs_sporadic, p_c_vs_sporadic

One-sided boundary LRT of CE vs sporadic (evidence of a familial effect at all).

chisq_ace_vs_ae, p_ace_vs_ae

One-sided boundary LRT of ACE vs AE (evidence C adds anything once A is estimated). A near-zero statistic with c2_ace at or near the 0 boundary indicates A and C are not jointly identifiable in this sample – report the simpler AE model (see package vignette / Leocadio-Miguel et al. companion analyses).

covs_tested, covs_kept, covs_dropped

Present only if covs was supplied and screen = TRUE: the full candidate list, the ones retained (p < prob_level), and the ones dropped.

covariate_lrt_p

Present only if covs was supplied and screen = TRUE: named numeric vector of each candidate covariate's individual LRT p-value against the full-covariate AE model.

Returns NULL if n < min_n.

Details

Unlike herit_vc(), which diagonalises a single variance component via eigen(), this function estimates one or two variance components jointly by direct numerical maximum likelihood over the full n x n covariance matrix (Cholesky-based GLS), since A and household indicator C do not in general share eigenvectors.

Model: Omega = sigma2_p * [h2*A + c2*C + (1 - h2 - c2)*I], jointly estimated for the ACE model via multi-start Nelder-Mead on a logit-reparametrised (h2, c2) (guaranteeing h2, c2 >= 0 and h2 + c2 < 1); AE and CE are 1-D special cases. Covariates (if any) enter the mean structure: E[y] = X %*% beta, with X including an intercept plus the (standardised) columns named in covs.

Covariate screening (when covs is given and screen = TRUE): matches SOLAR Eclipse's polygenic -screen, confirmed from source (lib/solar.tcl proc polygenic) rather than assumed. SOLAR's algorithm is a single pass, not iterative backward elimination: fit the AE model with all candidate covariates, then for each covariate individually, fit the AE model with just that one covariate removed and compute a likelihood-ratio test against the full-covariate model (chi-squared(1), not boundary-corrected, since a covariate coefficient is not a boundary parameter). Any covariate with p >= prob_level is dropped; all remaining covariates are kept. This happens once, against the AE model specifically (not separately for CE/ACE) – SOLAR's own script can in principle screen separately per model context, but using one shared covariate set for all four nested models here keeps the sporadic/AE/CE/ACE comparisons validly nested (a valid LRT requires the compared models to share the same mean structure).

LRTs for h2/c2: one-sided chi-squared(1) boundary correction, as in herit_vc(), since both null hypotheses (c2 = 0, or h2 = 0) sit on the boundary of the parameter space.

See also

build_grm(), build_household(), herit_vc(), herit_bivar()