Function map and cheat sheet
Source:vignettes/articles/function-map-cheatsheet.Rmd
function-map-cheatsheet.Rmd
library(hsquared)
#> hsquared is an experimental 0.9.0 release (not production; CRAN availability is tracked separately). Default fitting needs Julia + HSquared.jl; preview with hs_control(engine = "validate").
#> What you may report is listed on Can I fit and report this? -- not in validation_status(). Intervals are not coverage-calibrated.
#> See https://itchyshin.github.io/hsquared/articles/current-limits.html
# n = 4 is a syntax demo, not a number for a paper.
ped <- data.frame(
id = c("sire", "dam", "off1", "off2"),
sire = c(NA, NA, "sire", "sire"),
dam = c(NA, NA, "dam", "dam")
)
dat <- data.frame(
id = c("sire", "dam", "off1", "off2"),
sex = c("m", "f", "m", "f"),
weight = c(42, 38, 40, 37)
)One page for “which function do I call now?”. The preview below uses the same tiny pedigree as the README and Getting started. Every later example either runs without Julia or is marked as needing the engine.
The five steps
| Step | You are asking | Call |
|---|---|---|
| 1. Specify | Is my data in the right shape? Is my model expressible? |
hs_data(), data_status(),
animal(), formula_status()
|
| 2. Preview | What will the engine actually be sent? |
hs_control(engine = "validate"),
model_spec()
|
| 3. Fit | Estimate the variance components | hsquared() |
| 4. Check | Did it converge, and may I report it? |
fit_diagnostics(), summary(), Can I fit and report this?,
backend_info()
|
| 5. Extract | Get the quantities for the paper |
heritability(), variance_components(),
breeding_values(), fixef(),
ranef()
|
Steps 1, 2, and 4’s status calls need R only. Step 3 and the extractors need a fitted object, and therefore Julia.
1. Specify
# bundle phenotypes, pedigree, and any genomic/expression/environment metadata
bundle <- hs_data(phenotypes = dat, pedigree = ped)
summary(bundle)
#> <summary_hs_data>
#> components: phenotypes, pedigree
#> phenotype IDs: 4
#> ID overlap:
#> metric count
#> phenotype_ids 4
#> pedigree_ids 4
#> genotype_ids 0
#> expression_ids 0
#> phenotypes_without_pedigree 0
#> phenotypes_without_genotypes 4
#> genotypes_without_phenotypes 0
#> phenotypes_without_expression 4
#> expression_without_phenotypes 0
#> pedigree status:
#> metric count
#> pedigree_rows 4
#> pedigree_ids 4
#> phenotype_ids_with_pedigree 4
#> pedigree_only_ids 0
#> founders 2
#> nonfounders 2
#> known_sire_links 2
#> known_dam_links 2
#> missing_known_parent_ids 0
#> duplicate_pedigree_ids 0
#> self_parent_rows 0
#> same_known_parent_rows 0
data_status(bundle) # pedigree coverage, founders, parent links, marker map
#> <hs_data_status>
#> components: phenotypes, pedigree
#> ID overlap:
#> metric count
#> phenotype_ids 4
#> pedigree_ids 4
#> genotype_ids 0
#> expression_ids 0
#> phenotypes_without_pedigree 0
#> phenotypes_without_genotypes 4
#> genotypes_without_phenotypes 0
#> phenotypes_without_expression 4
#> expression_without_phenotypes 0
#> pedigree status:
#> metric count
#> pedigree_rows 4
#> pedigree_ids 4
#> phenotype_ids_with_pedigree 4
#> pedigree_only_ids 0
#> founders 2
#> nonfounders 2
#> known_sire_links 2
#> known_dam_links 2
#> missing_known_parent_ids 0
#> duplicate_pedigree_ids 0
#> self_parent_rows 0
#> same_known_parent_rows 0
#> expression status: not available
#> genotype status: not available
#> marker status: not available
#> annotation status: not available
#> environment status: not available
# what does the formula grammar parse, reserve, or refuse?
formula_status()
#> <hs_formula_status>
#> default: animal(); cbind()
#> opt-in:
#> permanent(); common_env(); maternal_genetic(); (1 | group);
#> animal(rr()); metafounder(); relmat(); precision(); genomic();
#> single_step()
#> reserved/planned: 14 reserved, 5 planned (error before fitting)
#> fences: $current_behavior
#> term
#> animal(1 | id, pedigree = ped)
#> animal(1 | id) with data = hs_data(..., pedigree = ped)
#> permanent(1 | id)
#> common_env(1 | group)
#> maternal_genetic(1 | dam)
#> (1 | group)
#> animal(rr(covariate, order = 2) | id, pedigree = ped)
#> maternal_env(1 | dam)
#> paternal_genetic(1 | sire, pedigree = ped)
#> paternal_env(1 | sire)
#> group(1 | genetic_group)
#> unknown_parent_group(1 | upg)
#> metafounder(1 | id, pedigree = ped, group = group, Gamma = Gamma)
#> inbreeding(1 | id)
#> cytoplasmic(1 | maternal_line)
#> imprinting(1 | id, pedigree = ped, parent = "maternal")
#> dominance(1 | id, pedigree = ped)
#> epistasis(1 | id, pedigree = ped)
#> relmat(1 | id, K = K)
#> precision(1 | id, Q = Q)
#> genomic(1 | id, Ginv = Ginv)
#> genomic(1 | id, markers = M)
#> single_step(1 | id, Hinv = Hinv)
#> single_step(1 | id, pedigree = ped, markers = M)
#> single_step(1 | id) with data = hs_data(..., pedigree = ped, genotypes = M)
#> single_step(1 | id, pedigree = ped, markers = M, group = group, Gamma = Gamma)
#> markers(M, model = "random")
#> marker_scan(M, map = marker_map)
#> qtl_scan(position, genotype_probs = probs)
#> cbind(trait1, trait2) ~ animal(1 | id, pedigree = ped)
#> animal(trait | id, pedigree = ped, cov = us())
#> animal(trait | id, pedigree = ped, cov = diag())
#> animal(trait | id, pedigree = ped, cov = lowrank(K = 2))
#> animal(trait | id, pedigree = ped, cov = fa(K = 2))
#> missing = miss_control(response = "include")
#> mi(x) with missing = miss_control(predictor = "model")
#> syntax_status fitting_status
#> parsed fitted (v0.1 default)
#> parsed fitted (v0.1 default)
#> parsed fitted (opt-in repeatability)
#> parsed fitted (opt-in common-environment)
#> parsed fitted (opt-in maternal)
#> parsed fitted (opt-in multi-effect)
#> parsed fitted (opt-in random-regression)
#> reserved not available
#> reserved not available
#> reserved not available
#> reserved not available
#> reserved not available
#> parsed fitted (opt-in supplied-Gamma metafounder)
#> reserved not available
#> reserved not available
#> reserved not available
#> reserved not available
#> reserved not available
#> parsed fitted (opt-in supplied relationship, experimental)
#> parsed fitted (opt-in supplied precision, experimental)
#> parsed fitted (opt-in genomic, covered validation-scale)
#> parsed fitted (opt-in genomic / SNP-BLUP)
#> parsed fitted (opt-in single-step)
#> parsed fitted (opt-in single-step construction)
#> parsed fitted (opt-in single-step bundle construction)
#> parsed fitted (opt-in supplied-Gamma H^Gamma)
#> reserved not available
#> reserved not available
#> reserved not available
#> parsed fitted (default route, covered multivariate)
#> planned not available
#> planned not available
#> planned not available
#> reserved not available
#> planned not available
#> planned not availableanimal(1 | id, pedigree = ped) is the additive genetic
term. With data = bundle you can write
animal(1 | id) and the bundle pedigree is used.
2. Preview: no Julia needed
hsquared(
weight ~ sex + animal(1 | id, pedigree = ped),
data = dat,
control = hs_control(engine = "validate")
)
#> Validated the v0.1 animal-model contract; no model was fitted (`engine = "validate"`). Julia fit target: `HSquared.fit_animal_model(y, X, Z, Ainv; method = :REML)`. The validated spec is returned invisibly as a named list; assign it to inspect the parsed contract (e.g. `spec$bridge$target`).
model_spec(
weight ~ sex + animal(1 | id, pedigree = ped),
data = dat
)
#> <hs_model_spec>
#> response: weight
#> family: gaussian(identity)
#> method: REML
#> observations: 4; fixed columns: 2; animal IDs: 4
#> animal term: animal(1 | id, pedigree = ped); group: id
#> random design: sparse Z with 4 nonzero entries
#> Ainv: build_in_julia
#> fitting: not run by `model_spec()`The preview reports reserved syntax before you start a Julia fit.
3. Fit: needs the engine
# needs Julia + HSquared.jl; not a silent fallback
fit <- hsquared(
weight ~ sex + animal(1 | id, pedigree = ped),
data = dat
)Opt-in experimental targets go through
hs_control(engine = "julia") with an
engine_control target. See Can I fit and report this? for which ones
you may report.
4. Check
summary(fit)
fit_diagnostics(fit) # convergence, iterations, at_boundary flag
backend_info() # which engine answered, and its version
validation_status() # developer evidence table, not the user listRead at_boundary before anything else. A variance pinned
at zero is a boundary result. Do not report it as an unconstrained
estimate.
5. Extract
| You want | Call |
|---|---|
| Heritability | heritability(fit) |
| Variance components | variance_components(fit) |
| Breeding values (EBVs/BLUPs) | breeding_values(fit) |
| Accuracy of those EBVs |
prediction_error_variance(fit),
reliability(fit), accuracy(fit)
|
| Fixed effects |
fixef(fit), coef(fit)
|
| Random effects | ranef(fit) |
| Fitted values and residuals |
fitted(fit), residuals(fit)
|
| Model comparison |
logLik(fit), AIC(fit) where supported;
nobs(fit)
|
| A plot |
plot(fit), autoplot(fit)
|
Route-specific extractors exist for the opt-in targets:
repeatability(), common_env_effects(),
maternal_effects(), genetic_covariance(),
marker_effects(), evolvability(), and
g_max(). ?g_matrix_geometry is the help topic
for that family; there is no g_matrix_geometry() function.
They inherit their route’s status.
Check fit_diagnostics() before comparing likelihoods.
Structured multi-effect and direct-maternal dense REML routes omit a
likelihood constant; AIC() refuses those results. A
matrix-free stochastic likelihood is available in diagnostics but not
through logLik() or AIC(). Compare fits only
when data, method, response family, and likelihood convention match. For
REML, the fixed-effect design must match as well.
R function and Julia function
R users start here. Some animal-model tools exist only in
HSquared.jl. The Julia pages map articles, not function
names; use this table for the function-level match.
| Task | R (hsquared) |
Julia (HSquared.jl) |
|---|---|---|
| Fit | hsquared() |
fit_animal_model(), fit_ai_reml()
|
| Summary of a fit |
summary(), print()
|
Julia only as fit_diagnostics(),
variance_components()
|
| h2, variance components |
heritability(), variance_components()
|
same names |
| h2 SE / interval |
heritability_standard_error(),
heritability_interval() (delta only) |
same names; also
heritability_interval(method = :profile)
|
| Other intervals | none | Julia only: variance_component_interval(),
bootstrap_variance_component_interval(),
genetic_correlation_interval()
|
Nested LRT / test V_A = 0
|
none (anova() stops; do not hand-build a
chi^2_1 LRT from logLik()) |
Julia only: nested_lrt() with
df = 1, boundary_df = 1 (50:50 chi^2_0 +
chi^2_1 mixture) |
| PEV, reliability, accuracy |
prediction_error_variance(),
reliability(), accuracy()
|
same names |
| A, A-inverse, inbreeding | none | Julia only: additive_relationship(),
pedigree_inverse(),
inbreeding_coefficients()
|
| Repeatability |
repeatability(),
repeatability_interval()
|
fit_repeatability_reml(),
repeatability_interval(); no
repeatability()
|
| Genomic relationship matrix | none (R takes Ginv or markers) |
Julia only: genomic_relationship_matrix(),
genomic_relationship_inverse()
|
| Marker scan |
gwas(); marker_scan() is reserved |
single_marker_scan(),
mixed_model_marker_scan(),
genome_wide_marker_scan()
|
| G-matrix geometry |
evolvability(), g_max(), and siblings;
help topic ?g_matrix_geometry
|
same extractor names |
Uncertainty: read this before using it
heritability_interval(fit) # experimental
variance_component_standard_errors(fit) # experimentalThese are asymptotic/delta-method quantities, surfaced when the
engine provides them and labelled experimental. They are not
coverage-calibrated. No route in this package currently carries
an interval-reporting permission. confint() and
vcov() are deliberately absent rather than shipped
uncalibrated.
Reserved names that do not fit yet
qtl_table(), gwas_table(),
eqtl_table(), lod_scores(), and
marker_variance_explained() are reserved output names for
future marker/QTL results. These names reserve vocabulary for later
work. marker_effects() is implemented on the opt-in
SNP-BLUP path.
Where to read next
- Getting started with hsquared: the guided introduction.
- Can I fit and report this?: reporting permissions per route.
- Fitting models: the longer fitting guide.